Goulnara N Arzhantseva

Victor S Guba
email: guba@uni-vologda.ac.ru

Martin Lustig
email: martin.lustig@univ.u-3mrs.fr

Jean-Philippe Préaux
email: preaux@cmi.univ-mrs.fr

J.-Ph Pr Éaux

Testing Cayley graph Densities (version longue)

Introduction

Let G be a group with a finite set of generators X of cardinality m. There is an associated Cayley graph C = C(G, X) (see §2), which has vertex set in bijection to G, and at every vertex there are precisely 2m adjacent edges. The combinatorial properties of the Cayley graph reflect the algebraic structure of the group G. In this paper we investigate the density of Cayley graphs, introduced first in [START_REF] Guba | On the properties of the Cayley graph of Richard Thompson's group F[END_REF]. This is a numerical parameter δ(C), defined below, which takes values between 0 and twice the number of group generators. It strongly depends on the isoperimetric properties of the Cayley graph and hence on those of G, which are often expressed in terms of the graph isoperimetric constant ι * (C) (see §2 for the definition). It is known that ι * (C) + δ(C) = 2m, see [START_REF] Guba | On the properties of the Cayley graph of Richard Thompson's group F[END_REF]. A group is amenable if and only if ι * (C) = 0, or, equivalently, δ(C) = 2m.

In order to estimate the density of a Cayley graph, one can compute densities of certain of its finite subgraphs. We propose a simple algorithm to construct an optimized subgraph (i.e. with a greater density) from any given finite subgraph of the Cayley graph. We apply the algorithm to amenable groups, non-amenable groups, and to groups for which it is not known whether they are amenable. More specifically, we investigate finitely generated free abelian groups, Baumslag-Solitar groups (amenable and non-amenable ones), the restricted wreath product Z Z, and Richard Thompson's group F . We analyze empirical data obtained by a C++ implementation of our algorithm.

We quote here only one of the numerical results obtained from our algorithm, which we find particularly interesting:

Sample result: There is a subset of cardinality 10169678 in Thompson's group F that has density 2.89577 with respect to the classical generating system of cardinality m = 2.

This work was supported by the Swiss National Science Foundation, No. PP002-68627.

Acknowledgements. Part of this work was conducted while the authors were visiting the Centre de Recerca Matemàtica in Barcelona. We would like to acknowledge the warm hospitality and support. Specifically we would like to thank José Burillo and Enric Ventura for their interest in our work, and for useful discussions.

Amenability and Følner families

Let G be a group generated by a finite set X. Let C = C(G, X) be the corresponding (left) Cayley graph. Recall that the set of vertices of C is G, and that the set of oriented edges is X ±1 × G. For any edge e = (x, g) the initial vertex is g, and the terminal vertex is xg. The inverse of the edge e, considered here separately from e, is the edge e -1 = (x -1 , xg). The label of e = (x, g) is defined to be the generator x ∈ X. The group G acts canonically on C from the right (by right multiplication of the vertices of C).

Throughout the paper we consider finite graphs A which are typically subgraphs of C. We always require that with any edge e also the inverse edge e -1 belongs to A.

The density of a finite graph A is defined by

δ(A) = v∈V (A) deg(v) #V (A) ,
where deg(v) denotes the number of oriented edges with initial vertex v, V (A) is the set of vertices of A, and #V (A) is the cardinality of V (A).

We define the density of the Cayley graph C = C(G, X) as supremum

δ(C) = sup A δ(A),
where A runs over all finite subgraphs of C.

Similarly, for any subgraph A of C one defines the isoperimetric constant

ι * (A) = #∂A #V (A) ,
where ∂A denotes the set of vertices of A that have adjacent edges in both, A and C -A. The infimum of the values of ι * (A), over all finite subsets of C, is called the isoperimetric constant of the graph C, and is denoted by ι * (C).

Lemma 1. For every finite subgraph A of C one has:

2m -(2m -1)ι * (A) ≤ δ(A) ≤ 2m -ι * (A).
Proof. Note that for any finite subset A ⊂ C the complement of ∂A in A consists entirely of vertices of degree 2m. Hence δ(A) = v∈V (A) deg (v) #V (A)

= v∈∂A deg(v) #V (A) + v∈V (A)-∂A 2m #V (A)
. But a vertex in ∂A has at least degree 1 and at most degree 2m -1 in A, which shows

#∂A #V (A) + 2m #V (A)-#∂A #V (A) ≤ δ(A) ≤ (2m -1) #∂A #V (A) + 2m #V (A)-#∂A #V (A)
and hence 2m -(2m -1)ι * (A) ≤ δ(A) ≤ 2m -ι * (A).

Isoperimetric properties of graphs play an important role in the study of amenable groups. There are many equivalent characterizations of amenability in the literature, see for example [START_REF] Ceccherini-Silberstein | Amenability and paradoxal decompositions for pseudogroups and for discrete metric spaces[END_REF] and the references given there. We use the following one.

Theorem 1. A finitely generated group G is amenable if and only if for some (or, equivalently, for any) finite generating set X the Cayley graph C = C(G, X) satisfies:

ι * (C) = 0

A family of finite subsets A n of C is called a Følner family (or a family of Følner sets if lim n→∞ ι * (A n ) = 0 .

In light of Lemma 1 this is equivalent to

lim n→∞ δ(A n ) = 2m .
Hence, the group G is amenable if and only if there exists a family of Følner sets A n ⊂ C.

For certain classes of groups there are well known Følner families. For example, if G is of polynomial growth, then one knows that with respect to any finite generating system X of G the set of balls B(n), which consists of all points in C of simplicial distance smaller or equal to n from the neutral element 1 ∈ G, is a Følner family [7,Proposition,Ch.VII.C.34]. Here we mean by simplicial distance the distance in the metric space obtained from C if one gives to every edge the length 1. Examples for groups of polynomial growth are free abelian groups and certain Baumslag-Solitar groups, which will be considered below.

However, many Baumslag-Solitar groups are of exponential growth, but some of them (not all !) are still amenable. The same is true for the wreath product Z Z, also considered below. In this case a Følner family exists in C, but the balls B(n) will not constitute such a family: There is a uniform upper bound strictly smaller than 2m to the density of every B(n). Of course, this last statement is true also if G is non-amenable.

Group presentations and normal forms

In order to compute in a finitely generated group G, one needs a normal from for the elements of G: For example, in Z × Z = a, b | aba -1 b -1 = 1 the element [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF][START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF] can be written as a 2 b, aba, baa, but also as aba -1 bab -1 a or a -69 ba 71 . It is an essential restriction on the class of groups G considered here that we require the existence of a uniquely determined normal form for the elements of G, and that this normal form can be recursively calculated. Notice that the generating set of G used in the normal form may well differ from the system X which is used to build the Cayley graph; in some cases this discrepancy is a rather convenient from a computational point of view.

3.1. Free abelian groups. The free abelian group of rank m is defined by the presentation:

x 1 , x 2 , . . . , x m | x i x j = x j x i for all i, j = 1, 2, . . . , m
A word in the canonical generators x 1 , x 2 , . . . , x m and their inverses is in normal form if and only if it is of the form

x p 1 1 x p 2 2 • • • x pm m
for some p 1 , p 2 , . . . , p m ∈ Z.

3.2.

Baumslag-Solitar groups. Let p, q ∈ Z \ {0}, and let BS(p, q) denote the Baumslag-Solitar group defined by the presentation:

a, b | ab p a -1 = b q
A word in a ±1 , b ±1 is in normal form whenever it is written in a reduced form (in the sense of HNN extensions [START_REF] Lyndon | Combinatorial group theory[END_REF]):

♠ b p 0 a q 1 b p 1 • • • a qn b pn
with n ≥ 0 and p 0 , p 1 , . . . , p n ∈ Z, q 1 , q 2 , . . . q n ∈ Z \ {0}, and such that for i = 1, 2, . . . , n -1 one has p i > 0, and p n ≥ 0. Furthermore, if q i > 0 then one has p i < |p|, and if q i < 0 then p i < |q|.

3.3.

The restricted wreath product Z Z. We define this group by the non-finite ♥ presentation:

a, x 0 , x 1 , . . . | x a i = x i+1 , x i x j = x j x i for all i, j > 0 A word in the generators or their inverses is in normal form whenever it is of the form

a n x p 1 i 1 x p 2 i 2 • • • x pn in , where n ∈ Z, 0 ≤ i 1 < i 2 < • • • < i n , and p 1 , p 2 , . . . , p n ∈ Z \ {0}.
In fact, this group can be generated by a and x 0 . We will refer to this as the canonical set of generators.

3.4.

Thompson's group F . Thompson's group F [START_REF] Cannon | Introductory notes on Richard Thompson's groups[END_REF] is the group of all piecewiselinear orientation preserving self-homeomorphisms of the unit interval such that (i) singular points are on dyadic numbers, and (ii) all slopes are integer powers of 2.

The group F admits the following infinite presentation:

x 0 , x 1 , x 2 , . . . | x j x i = x i x j+1 if i < j
It turns out that it has a finite presentation on two generators x 0 , x 1 . We will use this canonical set of generators in our computations. The generators x 0 and x 1 are given by the following functions, see also Figure 1.

x 0 (t) =    t/2 0 ≤ t ≤ 1/2 t -1/4 1/2 ≤ t ≤ 3/4 2t -1 3/4 ≤ t ≤ 1 x 1 (t) =        t 0 ≤ t ≤ 1/2 t/2 + 1/4 1/2 ≤ t ≤ 3/4 t -1/8 3/4 ≤ t ≤ 7/8 2t -1 7/8 ≤ t ≤ 1
We will consider two kinds of normal forms for the Thompson group F . The first one is given by words

x p 0 0 x p 1 1 • • • x pn n x -qn n • • • x -q 1 1 x -q 0 0
where n, p 0 , p 1 , . . . , p n , q 0 , q 1 , . . . , q n are non-negative integers such that (i) exactly one of p n or q n is non-zero, and (ii) if p k > 0 and q k > 0 for some 0 ≤ k < n, then p k+1 > 0 or q k+1 > 0. The left half x p 0 0 x p 1 1 • • • x pn n is called the positive part of the word and the right half

x -qn n • • • x -q 1 1
x -q 0 0 the negative part. A word is said to be positive (or negative) if its normal form only consists of its positive (or negative) part.

The second normal form is given by the so called reduced forest diagrams [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]. Recall that a binary forest is a finite sequence of binary trees, together with a pointer on one of the trees. The number of leaves in a binary forest is the sum of the numbers of leaves in its binary trees. A forest diagram is a pair of binary forests which have the same number of leaves. We speak of the bottom forest as well as of the top forest, see Figure 2.

Figure 2. A forest diagram with 8 leaves

A caret in a tree is a pair of leaves with the same parent vertex. A forest diagram is reduced if it has no opposite pairs of (bottom and top) carets (cf. [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]). For example, the diagram of Figure 2 is reduced.

One associates to an arbitrary reduced forest diagram an element of F with normal form

x p 0 0 x p 1 1 • • • x pn n x -qn n • • • x -q 1 1 x -q 0 0
as follows: -Enumerate top and bottom leaves, as well as top and bottom trees, from the left to the right, starting at 1. -The top (or bottom) forest gives the positive part (or negative part respectively) of the normal form. -The exponent of x i , for i > 0, equals to the maximal length of simple paths in the top forest starting at the i th top leaf and following the top-to-right direction (the exponent is 0 whenever such a leaf does not exist).

-The exponent of x -1 i , for i > 0, equals to the maximal length of simple paths in the bottom forest starting at the i th bottom leaf and following the bottom-to-right direction.

-The exponent of x 0 is n whenever the top pointer is on the (n + 1) st tree.

-The exponent of x -1 0 is n whenever the bottom pointer is on the (n + 1) st tree. For example, the reduced forest diagram of Figure 2 gives the element

x 0 x 2 1 x 2 4 x 6 x -1 7 x -1 6 x -3 2 .
Notice that adding a top and a bottom leaf on the right of a forest diagram does not change the corresponding element of F . Up to this trivial transformation, it turns out that each element of F can be represented by a unique reduced forest diagram (cf. [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]). This is the second normal form we are interested in.

For example, the generators x 0 , x 1 and their inverses are represented by the reduced forest diagrams given in Figure 3.

x 0 -1 x 1 -1 x 0 x 1 Figure 3.
Forest diagrams for the generators x 0 , x 1 and their inverses Below we use the following definitions: A binary forest is trivial if each of its subtrees consists of a single vertex only, and if the pointer is on the first of them. A forest diagram is negative (or positive) if its top (or bottom) tree is trivial; in such a case the normal form of the associated element of F is negative (or positive respectively). The height of a forest is the maximal height of one of its binary trees; it can take any value between 0 and the number of leaves minus 1. For instance, the top forest and the bottom forest of the forest diagram in Figure 2 have height 2 and 3 respectively.

Special subsets in Thompson's group F

For the first three classes of groups considered in this paper, free abelian groups, Baumslag-Solitar groups, and the wreath product Z Z, all of our numerical experiments are performed on balls B(n) of radius n in the Cayley graph C, centered around the vertex defined by the neutral element 1 ∈ G. For the fourth group, Thompson's group F , we will work with balls, but also with other kinds of sets, which we specify now. 4.1. Left positive balls. Let n > 0. The left positive ball of radius n, denoted by LP (n), is defined to be the maximal subgraph in the Cayley graph C = C(F, {x 0 , x 1 }) which contains only positive words

x p 0 0 x p 1 1 • • • x p k n with p 0 + p 1 + • • • + p k ≤ n
as vertices. Notice that a left positive ball LP (n) cannot be a tree: Indeed, our densification algorithm (see §5) deletes subtrees from some of the LP (n), but if n ≤ 5 the LP (n) are left unchanged (cf. §9.9). This shows in particular that those LP (n) are not trees. Thus none of the LP (n) is a tree, since obviously LP (n) is always a subgraph of LP (n + 1).

x 0 -1 x 1 -1 x 0 -1 x 2 -1 x 1 -1 x 2 -1 x 1 -2 x 2 -1 x 0 -1 x 0 -2 x 1 -1 e Figure 4. The negative forest N F (3)
negative forest with 3 leaves, N F (3), is given in Figure 4.

In the graphical representation of a negative forest we only draw the representative bottom forest. For example, in Figures 4 and5 below the generators x 0 and x 1 are given by simple and double arrows respectively.

Notice that even though N F (3) is a tree, N F (n) is not a tree for n ≥ 5 (see, for example, §9.10).

4.3.

Belk-Brown sets. The Belk-Brown set with n leaves and of height at most k, for n, k ∈ N and k < n, denoted by BB(n, k), is the maximal subgraph of the Cayley graph C which contains as vertices only elements that have a negative reduced forest diagram with n leaves and height at most k [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]. For instance, Figure 5 represents the Belk-Brown set BB [START_REF] Ceccherini-Silberstein | Amenability and paradoxal decompositions for pseudogroups and for discrete metric spaces[END_REF][START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF]. Note that BB(4, 1) is not a tree (it contains the loop

x -2 0 x -1 1 x 2 0 x -1 1 x -1 0 x 1 x 0 x 1 )
. Obviously one has BB(n, k) ⊂ BB(n + 1, k) and BB(n, k) ⊂ BB(n, k + 1). In particular BB(n, k) is not a tree whenever n ≥ 4 and k ≥ 1. Since the height of a forest is at most the number of its leaves minus 1, we have BB(n, n -1) = N F (n). This shows that N F (n) is not a tree whenever n ≥ 4, as already stated in the previous section.

The densification algorithm

In this section we describe the algorithm by which we can improve the density of a given finite graph, through passing over to a subgraph.

Given a finite subgraph A in the Cayley graph C of a finitely generated group G, the algorithm applies a finite sequence of reductions (given in detail below), and returns a new "densified" subgraph A. If the initial finite graph A is sufficiently

x 3 -1 x 1 -1 x 1 -1 x 2 -1 e x 0 -1 x 0 -3 x 0 -2 x 3 -1 x 0 -1 x 3 -1 x 0 -2 x 3 -1 x 3 -1 x 1 -1 x 0 -1 x 1 -1 x 0 -2 x 1 -1 x 0 -1 x 2 -1 x 0 -2 x 2 -1 x 0 -1 Figure 5. The Belk-Brown set BB(4, 1).
dense, then the returned graph A will have even higher density. Otherwise, for example if A is a cycle or a tree, it is possible that the algorithm will collapse A to a single vertex. As our applications concern all Cayley graphs built on 2-element generating systems, our algorithm is tuned to graphs with vertex degree uniformely bounded by 4. The necessary modifications for higher vertex degree, if needed, are fairly easy to device. We first need to introduce some terminology. Let K be a finite graph. A chain is a maximal simple path in K where all of its vertices, except for the endpoints, have degree 2 in K. The length of a chain is the number of its vertices of degree 2. A cycle is a simple loop in K where all vertices have degree 2 in K. A tripod is a subgraph of K which consists of 3 chains which have precisely 1 vertex v in common, and v is an endpoint in each of the three chains. The length of a tripod is the sum of the lengths of its three chains. A degenerated tripod is a subgraph of K which consists of two chains c and c , such that the two endpoints of c coincide with a vertex v, and precisely one of the endpoints of c also coincides with v. The length of a degenerated tripod is the sum of the lengths of the two chains c and c . We require that, in case of a tripod or a degenerated tripod, the vertex v has degree 3 in K; i.e. there is no other edge adjacent to v.

We will now define four types of elementary reductions. They correspond to the removal of chains, cycles, and tripods (degenerated or not) from K, whenever their length is large enough. Any such transformation will in most cases increases the density, see Lemma 2.

For a finite graph K with density δ(K) we define the following parameters: If (R3): Remove all chains of length greater than N c (K) from K. (R4): Remove from K some tripod of length greater than N t (K), or some degenerated tripod of length greater than N d (K). Repeat this procedure as often as possible, of course every time with respect to the new values of N t (K) and N d (K). Lemma 2. If K is a finite graph with density δ(K) > 2, then any of the above elementary transformations (R1), (R2), (R3) or (R4) transforms K into a subgraph K of strictly larger density of δ(K).

δ(K) = 2, let N c (K) = max(0, 2 δ(K)-2 ), N t (K) = max(0, 4 δ(K)-2 -1) and N d (K) = max(0, 2 δ(K)-2 -1). If δ(K) = 2, we set N c (K) = N t (K) = N d (K) = 0. ( R1
Proof. It suffices to check that the number of edges removed in any of the elementary transformations is strictly smaller than δ(K) times the number of vertices removed. This is trivially true for the transformations (R1) and (R2), since any tree has one more vertex than unoriented edges, and any cycle has equal number of vertices as unoriented edges. Since according to our conventions we have to count an edge and its inverse separately, this gives directly the desired inequality.

For the reduction (R3) we observe that any chain c of length n > N c (K) has precisely 2n + 2 edges (counting again an edge and its inverse separately) and n + 2 vertices, and while all of the edges are removed with c, only the n interior vertices of c are removed. Since δ(K) > 2 and

N c (K) = max(0, 2 δ(K)-2 ), one has n > 2 δ(K)-2 and hence δ(K) > 2 n + 2 = 2n+2 n ,
which precisely what we need. Any tripod of length n consists precisely of 2n + 6 edges and n + 4 vertices, where n + 1 of them (as well as all edges) will be removed. A degenerated tripod of length n consists precisely of 2n+4 edges and n+2 vertices, of which n+1 will be removed. The further calculation for (R4) is very similar to the above one for (R3) and thus left to the reader.

The algorithm proceeds as follows:

Given a finite graph K DO WHILE: the graph K is changing:

Apply successively reductions (R1), (R2), (R3), (R1), (R2), (R4) END RETURN \! the densified graph K Below we will call each successive applications of the reductions (R1), (R2), (R3), (R1), (R2), (R4) a round.

The algorithmic package

As mentioned in the introduction, our numerical results were obtained by means of computer calculations, executed by a program written in C++. In this section we give a brief description of the three parts I, II and III of our programmed algorithmic package, of its software routines and also of the assumptions and limitations involved. A fourth computational feature, concerning the linear interpolation of the numerical data obtained by parts I -III, and in particular calculating an interpolated limit density is performed using Matlab and is described in §8 below.

Part I of the algorithmic package consists of subprograms, one for each class of groups G considered here, that transform a given product of generators or their inverses into a word in normal form as introduced in §3.

Part II calculates, for a given parameter n, the finite graphs B(n), or, in case where the group in question is Thompson's F , the finite graphs LP (n), N F (n), or, for given k and n, the graph BB(n, k), as defined in §4.

Part III calculates, for any finite graph A (= the graph computed in part II) a densified subgraph A according to the algorithm presented in §5.

It is an important characteristics of our algorithmic package that it is organized in a strictly modular fashion where the different parts work independently from each other. This give the possibility to easily improve specific parts without having to change the rest. For example, new classes of groups can be investigated by adding new subprograms to part I without changing parts II or III, new families of sets for the known groups can be investigated without changing I or III, and the densification algorithm could be embellished and reapplied to the groups and set families already programmed without ever changing I or II.

6.1. What the program can do. The program works in a console mode. A contextual menu allows the user to choose any of the actions. The actions, besides saving a copy of the outputs into a text file and offering some further options, consists mainly in: -Choose one of the predefined groups; all further computations will concern this group.

-Perform direct computations, like writing an element or a product of elements in normal form.

-Construct one (or a sequence) of the predefined finite graphs as explained in §4, and compute their density (or alternatively their isoperimetric constant).

-Apply to such a finite graph (or sequence of finite graphs) the algorithm of §5.

The program provides some extra informations, like the density at each step, and further details concerning the application of the elementary reductions (R1)-(R4).

6.2.

What the program is made of. The program is written in standard C ++ and can be compiled either on Linux or W in32 platforms. It can be easily adapted to compilation on other platforms.

The program makes intensive use of the object-oriented abilities of C ++ . Groups, graphs, vertices are all objects (or class); all the main algorithms correspond to general functions which takes data as input (like a group), and returns data. The functions which construct balls in the Cayley graph, and implement the algorithm of §5 are general and can be applied to any implemented group or finite graph.

Elements are given by strings of characters. This allows more choices when one encodes an abstract group element. Usually, strings look like words on given canonical generators and inverses. So they are really close to their mathematical meaning. However, for example in the case of Thompson's group F , they don't represent words on the canonical generators, but encode normal forms. 6.3. Limitations. The graphs are constructed in the physical memory (RAM) of the computer, and their size is almost proportional to the number of vertices. The main limitation of our computation is obviously the size of the computed graphs. This is closely related to the complexity of the group: In groups of polynomial growth, our computations of B(n) can easily be implemented for n going up to hundreds or thousands, while in groups of exponential growth n goes hardly up to 20. On the other hand, the program can handle free groups of rank at most 128, free abelian groups of rank at most 128, and words in wreath products with i n at most 127. For the group F the program can handle normal forms x p 0 0 . . . x pn n x -qn n . . . x -q 0 0 with n up to 127.

Summary of experimental results

In this section we give an overview of our experimental results. For each of the classes of groups considered in §3 and for each type of the special subsets defined in §4 our presentation contains the following parts: (e) -comments. In the next section we give a graphical interpretation of our experimental results, and in Section 9 we present the numerical data obtained.

7.1. Comparative analysis I: Amenable vs non-amenable. (e) Balls in the Cayley graph of BS(1, -1) are isomorphic (as graphs) to balls in the Cayley graph of Z × Z. Hence one expects the same results as for Z × Z. This is indeed the case even if implementations of BS(1, -1) and of Z × Z are rather different: the former belongs to the class encoding Baumslag-Solitar groups, and the latter refers to the class encoding abelian groups. (e) This is the first example where the densification algorithm improves the (interpolated) limit density substantially. Of course, we know from the amenability of BS(1, 2) that there is some family of subgraphs which is Følner, and our calculations indicate that, even for large n, densification of balls B(n) is not sufficiently strong to build such a family.

Comparative analysis II:

The Thompson group F . Amenability of F is unknown, but one knows that F grows exponentially [START_REF] Cannon | Introductory notes on Richard Thompson's groups[END_REF], so that balls will certainly not give a Følner family. (e) The results are quite similar to the above case of the wreath product Z Z. An interesting point is that the densification algorithm performs at most three rounds for n < 15. However it suddenly takes 132 rounds to perform calculations for n = 15 and the density increases a lot. This allows to believe that there may well exist a subset of much higher density than given by our interpolation. (d) The interpolation of the limit density gives values 2.67 and 3.03 for the reduced negative forest and for the densified reduced negative forest respectively. (e) The interpolated values of the limit density are exceptionally close to the calculated ones: the norm of the residues is approx. 10 -5 . (d) The interpolation of limit densities gives 3.18, both before and after the densification. Thus the densification algorithm seems to be inefficient in this case. (e) Our interpolations do not agree with theoretical values: it has been announced that the limit of density of Belk-Brown sets tends to 3.5 [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF].

Notice that the best value of density of BB(n, k(n)) is obtained whenever k(n) increases, see Figures 40 and41. At the same time, the best value of densities of BB(n, k(n)) appear for k(n) = 3 (or for slowly growing k(n)).

For n fixed and k large enough all the BB(n, k) have the same density and the same number of vertices. It is an interesting question whether these finite graphs are isomorphic as subgraphs of the Cayley graph.

Graphics and interpolation

We give a graphical interpretation of our experimental results. The main numerical results are given in §9. Each subsection below concerns the density of a family of finite graphs in a given group, and its behavior under application of the densification algorithm. We study successively: (8.11) negative Forests N F (n) in F , (8.12) Belk-Brown's sets BB(n) in F .

(8.2) balls B(n) in Z × Z, (8.3) balls B(n) in Z × Z × Z, (8.4) balls B(n) in BS(1, -1), (8.5) balls B(n) in BS(1, 2), ( 8
Except for Thompson's group F , amenability (or not) of these groups is well known (see §2). In order to estimate the limit (or limit superior) of the density of the families of subgraphs considered, we apply a first order approximation to the numerical data obtained from our experiments. This interpolation allows us, in a certain sense, to extrapolate this limit of densities by a value called "interpolated limit density" of this family of subgraphs. Of course, the reader has to be aware that for groups with exponential growth this does only estimate a lower bound to the density of the Cayley graph, compare the discussion at the end of §2.

8.1. Method of interpolation. We consider the densities δ n of a sequence of finite graphs S n , for p + 1 ≤ n ≤ q. We approximate the δ n by a real-valued function f (n), specified below, defined on the domain {p + 1, p + 2, . . . , q}. We estimate the quality of the approximation first by constructing the vector in R q-p whose n th -component is the residue δ n -f (n), and later by considering its euclidian norm. We call this norm, the norm of residues. It is a non-negative number. Clearly, the smaller it is, the better the approximation will be, with the zero value for the norm of residues in the case of a perfect correspondence of f with the given values of δ n .

The approximating function f is set to be of type

f (n) = an + b n + x .
There are two main reasons: On one hand, these functions give the best experimental results. On the other hand, assume that the values δ n are well approximated by some rational function

f (n) = P (n) Q(n) .
Then one easily sees (observing that lim n→∞ δ n is neither zero nor infinite) that P and Q have to be of the same degree, and thus f (n) may be just as well approximated by a function of the above type an+b n+x . The key points for the interpolation procedure now are the following: First, the interpolation reduces to the consideration of a parameter x and a linear interpolation δ n (n+x) in order to obtain a and b. Here the value of x is chosen such that δ n (n+x) is best distributed close to a line. That is, the corresponding norm of residues in a linear interpolation is the smallest one among all possibilities for the value of x. (In the search for the best x, using M atlab, we consider only large enough n, and the values for x are only considered up to 10 -1 .)

Observe that the limit of f (n) for n → ∞ is a, and hence this is the value, called interpolated limit density, which we use as parameter to estimate the limit of δ n . The latter is, after all, the information we are mainly interested in. Aside: An interesting experimental result in the above described interpolation procedure is that the parameter a remains essentially unchanged whenever x is slightly modified. This stability with respect to perturbations seems interesting in light of the fact that a is related to the approximation of zero order: the line y = ax gives the asymptotic direction. 8.2. Free abelian group of rank 2. Numerical results are given in §9.1. We don't state here, although they were computed, large series of numerical data (up to n = 1000): they behave as expected. A first series of computation for δ(B(n)) presented below is going from n = 1 to 301, by laps of seize 10, and a second series is given for n = 1, 2, . . . , 15. This allows us to compare results for Z × Z with results for groups of exponential growth, and the "small scale" interpolation with the "large scale" interpolation as well as with the true behavior.

The "large scale" results for n = 1, 11, 21, . . . , 301 are given in Figure 6. We restrict ourselves to n = 11, 21, . . . , 301. The norm of residues for the density of balls B(n) is approx. 0.007, see Figure 7. The interpolation gives The interpolated limit density is equal to 4.

δ(B(n)) ≈ 4n -2 n + 0.5 .
The "small scale" results and interpolations, for n = 1, 2, . . . , 15, are given in Figure 9. For n = 2, 3, . . . , 15 the norm of residues of the interpolation of the density of the B(n) is approx. 0.018, see Figure 10. One obtains the interpolation δ(B(n)) ≈ 4n -1.4 n + 0.67 .

The interpolated limit density is equal to 4.

For n = 2, 3, . . . , 15 the densities of the densified B(n), multiplied by x + 0.67, give the line y = 4x -0.79. The norm of residues is approximately equal to 0.092, see Figure 11. The interpolation gives

δ(B(n)) = 4n 2 n 2 + n + 1/2 ∼ +∞ 4n n + 1 .
Our (small and large scale) interpolations are not exactly the same, but they are not so far off either. The key point is that in all three cases one obtains the limit density 4: The given subsets are a Følner family. The densities of the reduced balls is easy to deduce. The algorithm removes only four vertices and four edges (the 4 extremal points of the "square"):

δ(B(n)) = 4n 2 -4 n 2 + n -3/2 ∼ +∞ 4n n + 1
and the same conclusion holds.

8.3.

Free abelian group of rank 3. We proceed as above by comparing large scale and small scale interpolations with the true values. Numerical results are given in §9.2.

The "large scale" results for Z × Z × Z appear in Figure 12. Let us restrict to n = 11, 21, . . . , 171. The densities of balls B(n), rescaled via x + 0.6, are distributed along the line y = 6.x -5.4, see Figure 13. The norm of residues is approximately equal to 0.011. The interpolation gives

(n) in Z × Z × Z δ(B(n)) ≈ 6n -5.4 n + 0.6 ,
the interpolated limit density is equal to 6.

The norm of residues of the interpolation of the density of densified balls B(n) is approx. 0.032, see Figure 14. One obtains the interpolation

δ(B(n)) ≈ 6n + 4.3 n + 2.2 .
The interpolated limit density is equal to 6.

The "small scale" results for Z×Z×Z, for n = 1, 2, . . . , 15, are presented in Figure 15. We restrict ourselves to n = 3, 4, . . . , 15. y = 6.x -3.6, Figure 16. The norm of residues of the density of balls B(n) is less than 0.05, Figure 16. The interpolation gives:

δ(B(n)) ≈ 6n -3.6 n + 1
The interpolated limit density is equal to 6. As was expected, the computations give the same picture as in the case of Z × Z.

8.4. Baumslag-Solitar group BS(1, -1). Numerical results are given in §9.6.

The Cayley graphs of Z × Z and BS(1, -1) are not isomorphic. However, the underlying non-labelled graphs (given by the tiling of the plane by squares) are isomorphic. Thus our results, see Figure 18, are similar to those for Z × Z.

The interpolation of initial densities is given by: δ(B(n)) ≈ 4n -2 n + 0.5 The interpolated limit density is equal to 4. The interpolation of final densities is given by:

δ(B(n)) ≈ 4n + 1.61 n + 1.4
The interpolated limit density is equal to 4. The subsets B(n)) are to be Følner sets, as their density is minorized by that of the Følner family B(n). The interpolation is given by:

δ(B(n)) ≈ 3.22n -5.11 n -1.2
The interpolated limit density is equal to 3.22. 3.48n -11 n -2.9 , and the interpolated limit density is equal to 3.48.

8.6. Baumslag-Solitar group BS [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF][START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]. Graphics are given in Figure 22. Numerical results are found in §9.4. For n = 3, 4, . . . , 18, the density δ(B(n)), multiplied by (x -0.9), is distributed close to the line y = 2.64x -3.3, Figure 23. The norm of residues is approx. 0.066. The interpolation gives δ(B(n)) ≈ 2.64n -3.3 n -0.9 and the interpolated limit density is equal to 2.64. For n = 4, 5, . . . , 18, we multiply the density of the densified balls B(n) by x -1.9. This gives values that can be interpolated by the line y = 2.97x -6.29, see Figure 24. The norm of residues is approx. 0.084. Our approximation gives:

δ(B(n)) ≈
2.97n -6.29 n -1.9 The interpolated limit density is equal to 2.97. 8.7. Baumslag-Solitar group BS [START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF][START_REF] Cannon | Introductory notes on Richard Thompson's groups[END_REF]. Numerical results are given in §9.5, see Figure 25 for a graphical representation. For n = 3, 4, . . . , 14, the multiplied densities δ(B(n)) are distributed close to the line y = 2.44x -4.24, Figure 26. The norm of residues is less than 0.093. The interpolation gives δ(B(n)) ≈ 2.44n -4.24 n -1.6 , and thus the interpolated limit density is equal to 2.44. To interpolate the densities of the densified balls B(n), we restrict to values n = 3, 4, . . . , 14. The norm of residues is close to 0.26, see Figure 27. The interpolation is given by: δ(B(n)) ≈ 2.86n -5.08 n -1.5 The interpolated limit density is equal to 2.86. 8.8. Wreath product Z Z. Graphics are given in Figure 28. Numerical results are given in §9.7. The initial density, multiplied by x + 0.9, is distributed close to the line y = 2.43x + 0.55, see Figure 29. The norm of residues is approx. equal to 0.25.

The interpolation of densities is given by δ(B(n)) ≈ 2.43n + 0.55 n + 0.9 .

The interpolated limit density is equal to 2.43.

δ(B(n)) ≈ 3n -7 n -1.9
The interpolated limit density is equal to 3.

8.9. Balls B(n) in Thompson's group F . We now consider Thompson's group F , for which amenability is unknown. We first investigate the density of balls B(n) in F . Our results are given in Figure 31. Numerical results are given in §9. The estimation is given by: 2.8n -6.4 n -1.7 The interpolated limit density is equal to 2.8.

These results are clearly comparable with those of the amenable group Z Z.

8.10. Left-positive balls LP (n) in F . Numerical results are given in §9.9, graphical data are given in Figure 34. This gives the approximation:

δ(LP (n)) ≈ 2.22n -0.427 n + 0.343
The density of the densified sets LP (n), multiplied by n -0.7, is distributed close to the line y = 2.97x -6.25, see Figure 36. The norm of residues equals to 0.65841. This provides the approximation δ(LP (n)) ≈ 2.97n -6.25 x -0.7 .

The interpolated limit density is equal to 3.

8.11. Negative forests N F (n) in F . Numerical results are given in §9.10. A graphical interpretation is given in Figure 37. ).

n -1 n

For n = 5, 6, . . . , 14, the density of the densified negative forests N F (n), multiplied by n -1.1, is distributed along line y = 3.03x -6.28, Figure 39. The approximation is given by δ(N F (n)) = 3.03n -6.28 n -1.1 and the interpolated limit density is equal to 3.03. For k large enough all the densified BB(n, k) have the same density and the same number of vertices. Since BB(n, k) ⊂ BB(n, k + 1), one should expect that the sets obtained from our densification algorithm are the same. The best density we obtain is the density of BB(17, 3), which is equal to 2.89577.

We construct new sequences as follows. For each n we consider the maximum of densities of BB(n, k) as well as the maximum of densities of densified Belk-Brown sets BB(n, k). The sequences are given in Figure 42.

For n = 3, 4, . . . , 17 the best initial densities, multiplied with n+1.7, is distributed close to y = 3.18x -1.19, see Figure 43. The norm of residues is approx. 0.25. This

δ(n) = 3.18x -1.19 x + 1.7
The interpolated limit density is equal to 3.18.

For n = 5, . . . , 17, the best densities of the densified Belk-Brown sets BB(n, k), multiplied by n + 0.69, are close to the line y = 3.18x -2.92, see Figure 44. The norm of residues is less than 0.096. This gives as approximation:

δ(n) = 3.18x -2.92
x + 0.69 Thus, the interpolated limit density for the Belk-Brown sets BB(n, k) before densification and the one for the densified Belk-Brown sets BB(n, k) agree: they are both approximately equal to 3.18. In this section we will give some of the numerical data, obtained from our algorithmic package, in the form of tables. Each table corresponds to a fixed group and to a fixed family of subgraphs of the Cayley graph, as presented and discussed in the previous sections. For example the table given in §9.3 corresponds to the Baumslag-Solitar group BS [START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF][START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF] and to the family of balls B(n) in its Cayley graph. Each table is organised as follows:

A line in the table corresponds to a fixed choice of parameters for the family of subgraphs, thus specifying a particular subgraph. For example, the 5-th line of the table of §9.3 corresponds to the ball B( 4), for the group BS [START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF][START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]. The first colomn, labelled "Density", states the name of the subgraph considered. The second colomn, labelled "Before", states the density of this graph before applying the densification algorithm. The third colomn, labelled "#", states the number of rounds performed, i.e. the number of times the algorithm runs once through the sequence of steps specified at the end of §5 with at least one edge or vertex deleted (in which case BB (5,x) BB (6,x) BB (7,x) BB (11,x) BB( 12 it tries to repeat the maneuver, rising the value of # by one). The fourth colomn, labelled "After", states the density of the subgraph obtained as final result of our densification program. The fifth colomn, called "Increase", explicits the amount of density gained by the densification procedure, as well (in parathesis) the percentage this increase means with respect to the density before applying the program. Finally, the last colomn, labelled "Deleted vertices", states the number of vertices deleted in the densification procedure from the originally given subgraph. In parathesis it states the percentage this amounts to, with respect to the number of all vertices in the originally given graph.

Remark: If the densification algorithm is applied to a non-empty tree, then the number # of rounds the algorithm repeats the densification procedure will be equal to 1, and the final density must be equal to 0. This, however, is not specific for the case of non-empty trees: It will also happen, for example, if the given subgraph is a cycle, or any other graph obtained from gluing trees to a cycle.

Outlook

Throughout this section we use essentially the same notation as introduced in the beginning of the paper, and as used in the description of our densification algorithm ( §5): A always denotes a (not necessarily connected) finite graph, which we think of as "virtually" embedded as subgraph into an infinite ambient graph C. The latter is usually the Cayley graph C = C(G, X) of a finitely generated group G with respect to a generating system of finite cardinality m ∈ N, on which G acts on the right. In any case we always assume that C has a uniform bound 2m for the degree of any of its vertices. Note that the ambient graph is "virtual" in that, contrary to the subgraph A, it exists only as theoretical construct, and hence any finite piece A * of it (typically with A ⊂ A * ) has to be algorithmically constructed before it can be used in the algorithm.

For the purposes of this section it is easier to work with non-oriented edges. Thus every edge in this section corresponds to a pair of inversely oriented edges with same endpoints, in the notation of the earlier sections. Below we denote by V (A) the set of vertices of A and by E(A) the set of edges. By v(A) and e(A) we denote the cardinality of V (A) and E(A) respectively. We call e(A) the volume of A. The Euler characteristic of A is given by χ(A) = v(A) -e(A). For any vertex x ∈ V (A) the degree deg A (x) is the number of edge segments adjacent to x, which is consistent with the use of deg(A) in the earlier sections in view of the above transition from oriented to non-oriented edges.

A vertex x ∈ V (A) is called branch point if it has degree deg(x) ≥ 3.
For any subset X ⊂ A we denote by cl(X) the smallest subgraph of A that contains X. For any subgraph K of A we define the A-boundary ∂ A K of K to be the 0dimensional subgraph of K which consists of all vertices that bound simultaneously an edge from A and an edge from cl(A -K):

∂ A K = A ∩ cl(A -K) Also, int A K = K -∂ A K denotes the A-interior of K, which is in general not a subgraph. A subgraph K of A is called full, if it contains all edges of A that have both endpoints in K. The density of A is given by δ(A) = x∈V (A) deg A (x) v(A) = 2 e(A) v(A) .
We extend this notion in the obvious way to "graphs with some vertices missing" like the above set int A (K), for which one has δ(int ) . The comments and improvements proposed below concern the following three aspects: (A) the algorithmic determination of subgraphs of A with higher density, (B) the deterministic construction of larger graphs A * ⊂ C which contain A and have higher density, and (C) the non-deterministic construction of such A * .

A (K)) = 2 e(K) v(K)-v(∂ A K
A. Subgraphs with higher density. We observe that the improvements on the density by passing over to a subgraph A of A, as performed by the subroutines (R1) -(R4) of our algorithm presented in §5, are all based on the following principle: The computer checks for the existence of subgraphs K of A of a certain (fairly simple) type, and, if it finds any of them, it replaces A by A = A -int A K. The type of subgraphs K in question assures that the density increases strictly in this process. This is ensured by the topology of K, which needs to be of low density itself, and with small A-boundary. More precisely, one has:

Remark 1. For any subgraph K ⊂ A the complementary subgraph A = A -int A K satisfies δ(A) > δ(A) if and only if one has δ(A) > δ(int A K). Such a subgraph K is called density increasing.
Below we propose 4 further methods how to effectively find a density increasing subgraph K, in any given finite graph A:

(1) A first improvement of the algorithm used in our work can simply be obtained by embellishing the list of density increasing subgraphs K, which are integrated as fixed part of the algorithm without ever changing in the process. This is done by adding to the present list (i.e. trees, cycles, long chains, long tripods and long degenerated tripods, see §5) further subgraphs with low densities and small boundary. For example, any connected subgraph of A which is of class K(k, l, n), defined as set of all graphs K with v

(K) = k, |∂ A K| ≤ l and χ(K) ≥ n, is density increasing if k-n k-l < δ(A), by Remark 1. (2)
We devise a new subroutine, where the computer searches for the set K 0 of all vertices x ∈ V (A) with deg A (x) < δ(A) (or deg A (x) < δ(A) -C for some constant C > 0), and assembles them into "clusters", i.e. it builds iteratively full subgraphs K i which have a high percentage of low-valence vertices. The subgraphs K i are defined iteratively out of the connected components K j i-1 of K i-1 by adding vertices and edges from their neighborhood in order to create larger connected components, with the goal to decrease the total A-boundary of the union of the K j i-1 . At any given state K i the computer checks the cardinality of ∂ A K i , and stops the subroutine if this check shows that K i is density decreasing.

(3) A promising method to find interesting candidates for density decreasing subgraphs K ⊂ A seems to be the following: We consider a symmetric random walk on A where the starting measure on each vertex x is given by µ 0 (x) = 2m -deg A (x). We then let the random walk proceed for some integer time t, thus distributing the measure to give a value of µ t (x) on any x ∈ V (A) via the formula µ t (x) = µ t-1 (y) deg A (y) , where the sum is taken over all vertices y adjacent to x in A. For any h ≥ 0 we define the vertex sets V (h, t) = {x ∈ V (A) | µ t (x) ≥ h}, and the subgraphs K(h, t) as the full subgraphs of A with vertex set V (h, t). For any integer time t ≥ 0, if we let h decrease monotously from max{µ t (x) | x ∈ V (A)} to 0, the family K(h, t) defines a (finite) increasing nested sequence of subgraphs of A, which we propose as candidates for density decreasing subgraphs.

A variation of this approach would be to iterate the random walk until t is large enough so that the measure µ t (x) approximates a stable equilibrium µ ∞ (x), for all vertices x ∈ V (A). But this seems less interesting, as there is only one such limit ♠ distribution, and that is precisely given by 1 2m times the density function. Another, perhaps more promising variation comes from adding exterior measure sources or measure sinks, for example sinks for the high-density vertices, or sources for the low-density vertices, to force an equilibrium state to assemble the measure in the neighborhood of certain subgraphs considered as possible candidates for density decreasing subgraphs. (Recall that low-density does not imply density increasing, as one also needs that the A-boundary of K is small.) (4) More generally, improving the density of A by erasing interiors of subgraphs can also be viewed az improving the quality of A as an expander: We look for a "small" set Z of vertices (corresponding to ∂ A K in the above approaches) which cuts A into subgraphs K and A (with union A and intersection Z) that have rather different densities δ(A) > δ(K). If the difference of these densities is large with respect to the cardinality of Z, for example if δ(A) -δ(K) > 2m#Z v(K)-#Z , then δ(A) will be strictly bigger than δ(A).

There are also some interesting theoretical questions surrounding the algorithmic attempts to improve the density by erasing subgraphs: (i) Is there an algorithm to find the (possibly non-uniquely determined) subgraph A max of highest density among all subgraphs of A ?

As A is finite, the answer is of course "yes", but trying out all subgraphs is unfortunately not feasable in practise. Hence we rephrase the question as: (ii) What is the minimal complexity of any algorithm that derives A max from a given finite graph A. In particular, is there a polynomial-time algorithm ? (iii) Is there an algorithm for finding A max that uses only finitely many types of steps to pass from one intermediate subgraph A i to A i+1 ? Here a "step" consists of modifying a subgraph of A i of a given graph type into a new graph of given type. (iv) Is there always a sequence of nested subgraphs A i of increasing density and uniformly bounded volume difference e(A i ) -e(A i+1 ) connecting A to A max ? What is the minimal value for the volume difference bound needed to answer this question in the positive, in terms of the universal vertex degree bound 2m ?

An important fact the reader should note is the observation that certain "wrong" initial improvements on A (by erasing the interior of some density increasing subgraph K) can prevent the algorithm used in this paper, as well as any of the above proposed improvements ( 1) -( 4), from ever finding any of the really desired subgraphs A ⊂ A with density δ(A) close to δ(A max ). Indeed, it is not hard to find A i ⊂ A i+1 produced by the second program are automatically added to the list of test pairs K i ⊂ L i from L. The problem with this theoretically most promising approach is of course the hugh amount of memory needed to store the list L.

C. Non-deterministic methods for enlarging the test graph A. A different concept for finding high-density subgraphs of C comes from the observation that a random walk in a graph has the tendency to accumulate large amounts of measure (= "heat") in parts of the graph which are "heat preserving". There are several methods how to mimick random walks in more or less efficient ways on a computer: (1) For any non-negative function µ t : C → R which is equal to 0 outside a finite set A ⊂ C we define µ t+1 : C → R via µ t+1 (x) = µt(y) deg C (y) , where the sum is taken over all vertices y adjacent to x in C. We then define, for fixed t and increasing h, a decreasing family of finite level sets A t,h = {x ∈ C | µ t (x) > h} with empty intersection, which hence can be computed for any value of t and any value of h. Of course, the computation is very time consuming.

(2) An approximation of the previous method is the following: At any time t, one only distributes the weight µ t (y) among all of its neighbors x i (including possibly y itself) if µ t (y) is maximal or close to the maximum value of µ t (z) among all z ∈ A t = A t,0 .

The idea here is that vertices with small measure will have to be ignored anyway, as their totality grows too much like balls and will hence have low density, in general.

(3) We can exploit the fact that our graph C in question is not just any graph, but actually the Cayley graph of a group G, by denoting a (finite support) measure on C as element in the group ring RG. It is easy to see that convolution of (finite support) measures is nothing else than simply multiplying the corresponding elements in RG.

In particular, the classical nearest neighbor symmetric random walk on C(G, X) is directly given by the powers µ t for t → ∞, where µ = 1 2#X X∪X -1 x ∈ RG. (4) An interesting variation of the previous three approaches seems to be the following "discretization": One decides ahead of time on a finite integer scale, (say, from 0 to N ,) and rescale the heat function µ t at any time so that its maximal value on C equals N . Furthermore, for every vertex x the value µ t (x) is decreased to µ * t (x) = [µ t (x)], i.e. to the largest integer smaller or equal to µ t (x). This reduces on one hand the computational effort, and at the same time it cuts off the undesired very-low-heat vertices added by the pure random walk as described above in [START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF]. Of course, if one choses the scale too coarsely by picking N too small, we may get nowhere, by cutting off at every "rescaling second half-step" precisely what has been gained right before by the "neighbor-heat-distribution" in the first half-step, throughout any step of our discretized random walk procedure.

Figure 1 .

 1 Figure 1. The canonical generators x 0 and x 1 of F .

  ): Remove any subtree of K. (R2): Remove any cycle of K.

  (a) -known theoretical results for the group; (b) -best values of densities calculated by our program; (c) -analysis of the work of the densification algorithm; (d) -values of the interpolated limit density;

7. 1 . 1 .

 11 Free abelian groups Z × Z and Z × Z × Z. (a) Both groups are of polynomial growth, and hence balls B(n) are known to be a family of Følner sets. The slow growth allows an easy implementation of balls B(n) for large n (for hundreds or for thousands). Also, theoretical values of the density of balls and of densified balls are very easy to calculate. (b) The ball of radius 301 in Z × Z has density 3.98673 before and 3.98678 after the densification algorithm is applied. The denisities of the ball B(171) and of B(171) in Z × Z × Z are 5.94752 and 5.94812 respectively. (c) The densification algorithm does not change the initial density significantly: the increase of density is less than 1%. It deletes only 4 vertices in case of Z × Z.(d) The interpolated values of the limit density coincide with the theoretical values: 4 for Z × Z and 6 for Z × Z × Z respectively. Moreover, this is the case, both in small scale (for n = 1, . . . , 15) and in large scale (for n = 1, . . . , 300) calculations, which numerically confirms that balls constitute a Følner family. 7.1.2. Baumslag-Solitar group BS(1, -1).

  (a) The group is virtually abelian. Hence it is of polynomial growth and amenable, and balls B(n) are Følner sets. (b) The ball of radius 301 has density 3.98673. The densified ball of radius 301 has density 3.98678. (c) The densification algorithm removes only 4 vertices and the increase in density is negligible. (d) The interpolated limit density is equal to 4. In particular, one sees numerically that balls constitute a Følner family.

7. 1 . 3 .

 13 Baumslag-Solitar group BS[START_REF] Fordham | Minimal Length Elements of Thompson's Group F[END_REF][START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF]. (a) This is an amenable group of exponential growth. Hence balls do not form a Følner family. (b) The ball of radius 19 has density 3.14771. The densified ball of radius 19 has density 3.42439.(c) The densification algorithm yields a 9% increase in density and slightly reduces the size of balls: approx. 40% of vertices are removed. (d) Our interpolation gives limit density 3.22 for balls and 3.48 for the densified balls. They are quite close to the optimal value 4.

7. 1 . 4 .

 14 Baumslag-Solitar group BS[START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF][START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF].(a) This group contains F 2 × Z as a subgroup (of index 2) and hence it is not amenable. (b) The ball of radius 18 has density 2.58585. The densified ball of radius 18 has density 2.928. (c) The densification algorithm induces a 14% increase in density and remove approx. 40% of vertices. (d) The interpolated values of the limit density of balls and of densified balls are 2.64 and 2.97, respectively. (e) An interesting point is that the algorithm runs only through one round.7.1.5. Baumslag-Solitar group BS[START_REF] Belk | Forest diagrams for elements of Thompson's group F[END_REF][START_REF] Cannon | Introductory notes on Richard Thompson's groups[END_REF].

  (a) This group is not amenable. Hence the density of balls can not be close to the optimal value 4. (b) The ball of radius 14 has density 2.40677. The densified ball of radius 14 has density 2.79782. (c) The densification algorithm generates a 16% increase in density and removes relatively many vertices: approx. 70%. (d) The interpolated values of the limit density of balls and of densified balls are 2.44 and 2.86, respectively. 7.1.6. Wreath product Z Z. (a) This is an amenable group of exponential growth. Balls B(n) are not a Følner family. (b) The ball of radius 16 has density 2.32838. The densified ball of radius 16 has density 2.90938. (c) The densification algorithm produces a 25% increase in density and removes approx. 2/3 of vertices. Thus the algorithm is quite efficient in this case. (d) The interpolated values of the limit densities of balls and of densified balls are 2.43 and 3 respectively.

  7.2.1. Balls B(n) in F . (b) The ball of radius 15 has density 2.14905. The densified ball of radius 15 has density 2.7183. (c) The densification algorithm induces a 25% increase in density. It removes more than 80% of vertices. (d) The interpolated limit densities of balls and of the densified balls are 2.23 and 2.8 respectively.

7. 2 . 2 .

 22 Left positive balls LP (n) in F . (b) The left positive ball LP (19) has density 2.15988. The densified left positive ball LP (19) is of density 2.74349. (c) The densification algorithm yields an increase in density of approx. 27%. This is one of the best values obtained. The densified left positive balls are particulary small: up to 90% of the vertices are removed by the densification algorithm. (d) The interpolated limit densities of left positive balls and of densified left positive balls are 2.22 and 2.97 respectively. (e) The densification algorithm appears to be most efficient in case of these particular graphs. 7.2.3. Negative forests N F (n) for F . (b) The reduced negative forest N F (14) has density 2.47619. The densification of this subgraph gives density 2.79448. (c) The densification algorithm is rather inefficient: it gives a 13% increase in density and removes less than 60% of the vertices.

7. 2 . 4 .

 24 Belk-Brown sets BB(n) for F . (a) Since BB(n, n -1) = N F (n) both implementations have the same behavior, even though distinct routines and functions are used. (b) The Belk-Brown set BB(17, 3) has density 2.82642. The densified set BB(17, 2) is of density 2.89577. This is comparable to the case of Z Z. (c) The densification algorithm increases density by less than 13% and removes approx. 60% of vertices.

Figure 6 .Figure 7 .

 67 Figure 6. Large scale density of balls B(n) in Z × Z

Figure 8 .Figure 9 .

 89 Figure 8. Large scale interpolation of densities of densified balls B(n) in Z × Z

Figure 10 .Figure 11 .

 1011 Figure 10. Small scale interpolation of densities of balls B(n) in Z × Z

Figure 12 .

 12 Figure 12. Large scale density of balls B(n) in Z × Z × Z

Figure 14 .Figure 17 .

 1417 Figure 14. Interpolation of large scale density of densified balls B(n) in Z × Z × Z

Figure 18 .

 18 Figure 18. Density of balls B(n) for BS(1, -1)

8. 5 .

 5 Figure 19. Density of balls B(n) in BS(1, 2)

Figure 20 .

 20 Figure 20. Interpolation of densities of balls B(n) in BS(1, 2)

Figure 22 .

 22 Figure 22. Density of balls B(n) of BS(2, 2)

Figure 24 .Figure 25 .

 2425 Figure 24. Interpolation of densities of densified balls B(n) in BS(2, 2)

Figure 26 .Figure 27 .

 2627 Figure 26. Interpolation of densities of balls B(n) in BS(2, 3)

8 .

 8 In order to interpolate the densities δ(B(n)) we only consider n = 3, 4, . . . , 15. Multiplication of the density of B(n) by n + 2.3 gives values close to the line y = 2.23x + 3.65, see Figure 32. The norm of residues is approx. 0.123.

Figure 28 .

 28 Figure 28. Density of balls B(n) in Z Z

Figure 30 .Figure 31 .

 3031 Figure 30. Interpolation of densities of densified balls B(n) in Z Z

Figure 32 .

 32 Figure 32. Interpolation of densities of balls B(n) in F

Figure 34 .

 34 Figure 34. Density of left-positive balls LP (n) in F

Figure 36 .Figure 37 .

 3637 Figure 36. Interpolation of the densities of the densified sets LP (n) in F

Figure 38 .

 38 Figure 38. Interpolation of densities of N F (n) in F

Figure 40 .

 40 Figure 40. Densities of Belk-Brown sets BB(n, k) in F

Figure 41 .

 41 Figure 41. Final densities of Belk-Brown sets BB(n, k) in F

Figure 42 .

 42 Figure 42. Best densities of Belk-Brown sets

Figure 44 .

 44 Figure 44. Interpolation of best densities of the densified Belk-Brown sets BB(n, k)

  4.2. Negative forests.A negative forest with n leaves, n ∈ N, denoted by N F (n), is defined to be the maximal subgraph of the Cayley graph which contains only vertices that are given by group elements which are represented by a negative reduced forest diagram with at most n leaves. Obviously one has N F (n) ⊂ N F (n + 1). The

examples of graphs (for example built on two disjoint graphs connected by adding a long chain) which answer the following question in the negative:

Is any subgraph A of A with δ(A) ≥ δ(A), such that A does not contain a subgraph of strictly larger density than δ(A), necessarily equal to A max ? B. Deterministic methods do increase A to a larger graph A * with higher density. We first notice that for any finite subgraph A of the Cayley graph C of G, and for any g ∈ G with sufficiently large translation length in C, the subgraph Ag of C is a disjoint isomorphic copy of A, and hence their union has the same density as A. On the other hand, if one finds an element g ∈ G such that A and Ag intersect in a single vertex, then the density of the union A ∪ Ag is strictly larger than that of A. Of course, as the special case of a subtrees A shows, there are rather strict limits to this method in its crude form, but nevertheless it gives the right idea why the following is promising.

Since A is finite, the subset G A ⊂ G defined by

is also finite, so that at least in principle one can calculate, for all subsets B ⊂ G A , the density of the union

It seems quite realistic that among the A B one finds new graphs with substantially higher density than A, and that an iteration of this procedure leads to a very promising family of density test graphs for G.

On the other hand, the calculation of the unions A B is tedious and requires much computing time. Hence the following suggestion may prove to be helpful: Denote Ag by A , and let K = A∩A be the intersection subgraph. Then δ(A∪A ) is calculated by the formula

and hence δ(A ∪ A ) > δ(A) if and only if δ(A) > δ(K).

We may thus start out with a large g = g q = x q . . . x 1 ∈ G, so that Ag is disjoint from A and then pass successively to Ag q-1 , to Ag q-2 , etc, for g k = x k . . . x 1 , until A and Ag k meet. As small graphs have (a forteriori) small density, the first non-empty intersection graphs K = A ∩ Ag i seem to be interesting candidates for the above procedure.

A very different deterministic approach to construct families A i of increasing volume e(A i ) and increasing density δ(A i ) consists of systematic "local" improvements implemented as follows:

A first computer program compiles a complete list L = L n , for some integer n ≥ 1, of all pairs of subgraphs K i ⊂ L i contained in the ball B n (1) in C of radius n around the trivial element 1 ∈ G, which satisfy δ(L i ) > δ(K i ). A second program then verifies, for any x ∈ A i , whether B n (x) ∩ A i = K i x, and if so, replaces the subgraph K i x of A i by L i x to obtain the new graph A i+1 .

Of course, if one can increase the index n of the list L n in the first computer program, then procedure performed by the second program will lead to better values. In principle one can also imagine an interactive procedure, where all pairs of graphs


Jacek Bieroń
email: bieron@uj.edu.pl

Paul Indelicato
email: paul.indelicato@spectro.jussieu.fr

Per Jönsson
email: per.jonsson@ts.mah.se

Multiconfiguration Dirac-Hartree-Fock calculations of transition rates and lifetimes of the eight lowest excited levels of radium

The multiconfiguration Dirac-Hartree-Fock (MCDHF) model has been employed to calculate the transition rates between the nine lowest levels of radium. The dominant rates were then used to evaluate the radiative lifetimes. The decay of the metastable 7s7p 3 P0 state through 2-photon E1M1 and hyperfine induced channels is also studied.

Introduction

Recent advances in trapping and spectroscopy of free, neutral atoms make it possible to extend the search for time reversal violation effects into the domain of radioactive elements [1]. In the last decade several heavy atoms were considered as candidates for experimental searches [2]. There are at least two ongoing atomic trap experiments (in Kernfysisch Versneller Instituut in Netherlands [3,4] and in Argonne National Lab in the U.S. [5,6]) whose aim is to detect the electric dipole moment of radium. The advantage of radium lies in octupole deformations of nuclei in several isotopes [7], simple electronic structure (ground state configuration [Kr]4d 10 4f 14 5s 2 5p 6 5d 10 6s 2 6p 6 7s 2 yields the closed-shell singlet state 1 S 0 ) as well as in coincidental proximity of two atomic states of opposite parity, 7s7p 3 P 1 and 7s6d 3 D 2 , which are separated by a very small energy interval 5 cm -1 .

The data on atomic spectrum of radium compiled in the tables of Moore [START_REF] Moore | Atomic Energy Levels, NSRDS-NBS No[END_REF] came from the experimental investigations of Rasmussen [START_REF] Rasmussen | [END_REF], with subsequent revisions by Russell [10]; both go back to the 1930s. These data cover only 69 classified lines. The isotope shifts and hyperfine structures of radium were measured by the group of Wendt [11,12,13,14,15,16] in the 1980s. They have studied both atomic Ra I and ionic Ra II spectra, and obtained the isotope shifts, magnetic dipole hyperfine constants A and electric quadrupole constants B of the 7s7p 1 P 1 , 7s7p 3 P 1 , 7s7p 3 P 2 , and 7s7d 3 D 3 levels of neutral radium, as well as the magnetic dipole hyperfine constants A of the 7s 2 S 1/2 , 7p 2 P 1/2 , and 7p 2 P 3/2 levels of Ra + ion, together with the electric quadrupole constant B of the 7p 2 P 3/2 state of Ra + . Hyperfine structures of singly-ionised radium have been also the subject of several theoretical papers [17,18,19,20,21,22].

The excitation energies of several states of neutral and singly-ionised radium (and barium) were later calculated in the framework the relativistic coupled-cluster theory by the group of Kaldor [23,24]. More recently there have been three papers from our group in which we calculated the lifetimes of the 7s7p 1 P 1 and 7s6d 3 D 2 states [25]; the hyperfine structure constants of all levels belonging to the two lowest excited-state configurations [26]; and the electric field gradients generated by the electronic cloud in the 7s7p 1 P 1 , 7s7p 3 P 1 , and 7s7p 3 P 2 states, which in turn (combined with the measured values of the electric quadrupole constants B) yielded the nuclear electric quadrupole moment of radium-223 isotope [27]. The excitation energies and lifetimes of several states of radium (and barium) were calculated recently in the framework of the combined configuration-interaction and many-body perturbation theory by Dzuba and Ginges [28].

In the present paper we calculated the transition probabilities between the states arising from the three lowest configurations of radium: 7s 2 1 S 0 , 7s7p 3 P 0,1,2 , 7s7p 1 P 1 , 7s6d 3 D 1,2,3 , and 7s6d 1 D 2 , as well as the lifetimes of these states. The purpose of the present paper is fourfold. Firstly, we intended to extend the transition rate calculations on all levels arising from the three lowest configurations (i.e. 7s 2 , 7s7p, and 7s6d). We included (1) transitions which contribute appreciably to lifetimes; (2) transitions involving 7s6d 3 D 2 and 7s7p 3 P 1 levels because they are of interest in EDM experiments (the mixing induced by T-odd interactions is strongest between these two levels); (3) transitions which may be important for trapping; (4) transitions which are stronger than 1/s (a somewhat arbitrary threshold). Secondly, new comparison of our results became available when the two most recent papers [28,5] appeared in print. These in turn permitted further tests of the newly developed [29] parallel version of the GRASP package [30] and calibration the theoretical model for the calculations of the spectroscopic properties of radium. Finally, we present a summary of the available theoretical data on the spectroscopic properties of eight lowest excited levels of radium, with the hope that they may be of help for the experimental groups, that are currently in the process of setting up the atomic traps for the search of permanent electric dipole moments. The Argonne group has already measured [5] the frequency and the rate of the 3 P 1 -1 S 0 transition in the 225 88 Ra isotope, and determined the lifetime of the 3 P 1 level. The lifetime of this level have been previously calculated by Hafner and Schwarz [31] and by Bruneau [32].

Theory

The modified version [29] of the GRASP implementation [30] of the multiconfiguration Dirac-Hartree-Fock method [33] was used in the present paper. The starting point is the Dirac-Coulomb Hamiltonian

H DC = i cα i • p i + (β i -1)c 2 + V N i + i>j 1/r ij , (1)
where V N is the monopole part of the electron-nucleus Coulomb interaction. The wavefunction for a particular atomic state (Ψ ) is obtained as the self-consistent solution of the Dirac-Fock equation [30] in a basis of symmetry adapted configuration state functions (Φ)

Ψ (Γ P JM ) = N CF i c i Φ(γ i P JM ). ( 2
)
The basis N CF was systematically enlarged [34,35] to yield increasingly accurate approximations to the exact wavefunction. All calculations were done with the nucleus modeled as a variable-density sphere, where a two-parameter Fermi function [36] was employed to approximate the charge distribution. The Breit and QED corrections were estimated with the step-wise procedure described in [34]. They were applied only to two transition rates, as explained in section 4.3 below.

Method

The wavefunctions were obtained with the active space method in which configuration state functions of a particular parity and symmetry are generated by substitutions from a reference configuration to an active set of orbitals. The active set and the multiconfiguration expansions are increased systematically. The whole process is governed by convergence of the expectation values. The calculations were divided into two stages. Each stage was further divided into several consecutive steps. In the first stage, the spectroscopic and virtual spinorbitals were generated in relatively small multiconfiguration expansions. The spectroscopic orbitals required to form a reference wavefunction were obtained with a minimal configuration expansion, with full relaxation. Then virtual orbitals were generated in five consecutive steps. At each step the virtual set has been extended by one layer of virtual orbitals. A layer is defined as a set of virtual orbitals with different angular symmetries. In the present paper five layers of virtual orbitals of each of the s, p, d, f, g, h symmetries were generated. At each step the configuration expansions were limited to single and double substitutions from valence shells to all new orbitals and to all previously generated virtual layers. These were augmented by small subsets of dominant single and double substitutions from core and valence shells, with further restriction, that at most one electron may be promoted from core shells (which means, that in the case of a double substitution the second electron must be promoted from a valence shell). All configurations from earlier steps were retained, with all previously generated orbitals fixed, and all new orbitals made orthogonal to others of the same symmetry. The initial shapes of radial orbitals were obtained in Thomas-Fermi potential, and then driven to convergence with the self-consistency threshold set to 10 -8 . All radial orbitals were separately optimized for each of the nine atomic states of interest. The Optimal Level form of the variational expression [36] was applied in all variational calculations. In the second stage, the configuration-interaction calculations (i.e., with no changes to the radial wavefunctions) were performed, with multiconfiguration expansions tailored in such a way, as to capture the dominant electron correlation contributions to the expectation values. The valence and core-valence effects constitute the dominant electron correlation contributions in the oscillator strength calculations [START_REF] Fischer | Computational Atomic Structure. An MCHF Approach[END_REF], therefore all single and double substitutions were allowed from several core shells and both valence shells (i.e., 7s 2 , 7s7p, or 7s6d, depending on the state) to all virtual shells, with the same restriction as above, i.e. that at most one electron may be promoted from core shells. The virtual set was systematically increased from one to five layers, until the convergence of transition rates was obtained. In a similar manner, several core shells were systematically opened for electron substitutions -from the outermost 6p to 5s5p5d6s6p shells. The effects of substitutions from 4s4p4d4f shells were neglected. We estimated them separately for three states and discovered that they change the calculated values of transition rates by no more than a fraction of a percent. The transition rates were calculated with the biorthonormal technique [START_REF] Malmqvist | [END_REF]39], which permits the application of standard Racah algebra, while retaining the advantage of wavefunctions separately optimized for each state. Experimental values of transition energies from Moore's tables [START_REF] Moore | Atomic Energy Levels, NSRDS-NBS No[END_REF] were used in calculations of transition rates.

Results

4.1 The metastable 7s7p 3 P 0 state In principle there are three possible decay channels of the 7s7p 3 P 0 state. It can decay to (the only lower lying) ground state through (1) a blackbody radiation induced decay, (2) a 2-photon E1M1 transition, or through (3) a hyperfine induced transition. The first is beyond the scope of the present paper, since it depends on the ambient temperature. Of the other two, the former can be estimated through an order-of-magnitude comparison with the E1M1 1s2p 3 P 0 -1s 2 1 S 0 two-photon transition in helium-like heavy ions. We started from a recent evaluation [40], which gives 3.14 × 10 9 s -1 for Ra. In order to obtain a good dependence on transition energy and radial matrix elements, we evaluated the E1M1 matrix elements (following Eq. (2) of Ref. [40]) for He-like Ra (using only 1s2p 3 P 1 as an intermediate state) and for neutral Ra (using only the 7s7p 3 P 1 ). Then we integrated over photon energies with the five-point Gauss-Legendre formula. Finally, we evaluated the ratio of these two values to scale the He-like E1M1 rate. This gives an order of magnitude estimate of 9.6 × 10 -3 s -1 . Such lifetime (≈ 100 s) of the 7s7p 3 P 0 state would be comparable to the lifetimes of nuclei of several radium isotopes (it would be significantly shorter only in comparison with the nuclei of the most stable radium isotopes spanning the mass range 223-229). This would be the case of spin-zero isotopes of radium.

For the isotopes of radium with a nonzero value of nuclear spin, the hyperfine-induced transition must also be considered. We estimated the 7s7p 3 P 0 -7s 2 1 S 0 transition rate with a simple three-state model, in which the wavefunctions of the hyperfine components of the upper 3 P 0 state are described by a symmetry-adapted configuration-state-function expansion of the form

|7s7p 3 P 0 IF HF S = c 0 |7s7p 3 P 0 IF + c 1 |7s7p 1 P 1 IF + c 2 |7s7p 3 P 1 IF (3)
running over the appropriate hyperfine components IF of the 1 P 1 and 3 P 1 states. The hyperfineinduced transition rate (in s -1 ) may be approximately expressed as

A(7s7p 3 P o 0 → 7s 2 1 S 0 ) = 2.02613 × 10 18 3λ 3 c 1 7s7p 1 P o 1 Q (1) 1
7s 2 1 S 0 + c 2 7s7p 3 P o 1 Q (1) 1 7s 2 1 S 0 2 (4)
where 7s7p 1 P o 1 Q

(1) 1

7s 2 1 S 0 and 7s7p 3 P o 1 Q

(1) 1

7s 2 1 S 0 are reduced matrix elements for the electric dipole operator, and λ is the transition wavelength (in Å). The coefficients c 1 and c 2 are related to the off-diagonal magnetic dipole constants A HF S M1 ( 1 P 1 , 3 P 0 ) and A HF S M1 ( 3 P 1 , 3 P 0 )

c 1 = I(I + 1) A HF S M1 ( 1 P 1 , 3 P 0 ) ∆E( 1 P 1 -3 P 0 ) , c 2 = I(I + 1) A HF S M1 ( 3 P 1 , 3 P 0 ) ∆E( 3 P 1 -3 P 0 ) (5)
(see [41] or [42] for full derivation). The sum in the expansion (3) should run over all excited states, but the sum in Eq. ( 4) is usually dominated by those states, for which the transition rates to the ground state are large and at the same time the energy denominators in Eq. ( 5) are small. In case of the 7s7p 3 P 0 state of radium the 7s7p 3 P 1 and 7s7p 1 P 1 states dominate. The contribution of the reduced matrix element 7s7p 3 P o 1 Q

(1) 1 7s 2 1 S 0 is 3.7 times larger than that of the 7s7p 1 P o 1 Q

(1) 1

7s 2 1 S 0 matrix element. The contributions of other states are much smaller, due to the presence of the energy denominators in the coefficients (5) and to the fact that the corresponding off-diagonal matrix elements are many orders of magnitude smaller (actually they are exactly zero, unless non-orthogonality between fully relaxed wavefunctions and correlation effects are not neglected), and they can be safely ignored at present level of the overall accuracy. The rate of the hyperfine-induced transition is isotope-dependent, i.e. it depends on the nuclear spin and on the nuclear magnetic moment. As in our previous paper [27], the 223 88 Ra isotope was chosen to set the nuclear parameters (the transition rate A( 3 P 0 -1 S 0 ) may be readily recalculated for other radium isotopes, for which nuclear spins and magnetic moments are known). The nuclear spin of 223 88 Ra is I = 3/2, and the nuclear magnetic dipole moment µ = 0.2705 (19)µ N was taken from the paper of Arnold et al [13]. The electric dipole transition rates A( 1 P 1 -1 S 0 ) and A( 3 P 1 -1 S 0 ) from table 1 (in Babushkin gauge) were used. The off-diagonal magnetic dipole constants A HF S M1 ( 1 P 1 , 3 P 0 ) = 540 MHz, and A HF S M1 ( 3 P 1 , 3 P 0 ) = 1172 MHz, were evaluated with the use of the same wavefunctions, and in the same approximation, as described in section (3) above. Together they yield hyperfine-induced transition rate A( 3 P 0 -1 S 0 ) = 0.0210 s -1 in case of constructive interference, and A( 3 P 0 -1 S 0 ) = 0.0070 s -1 in case of destructive interference.

The above evaluations of the E1M1 and hyperfine-induced rates were performed independently. As it turned out, the two contributions are of the same size, therefore they have to be treated simultaneously. To this end we employed the effective Hamiltonian method [43]. It requires also an evaluation of the A HF S M1 ( 3 P 1 , 1 P 1 ) matrix element (851 MHz). The method is valid beyond the limits of the perturbation method exposed above, in particular it does not require that the energy separation between levels is large compared to level widths. It yields A 3 P 0 -1 S 0 = 0.02935 s -1 (if the two-photon transition is switched off the effective Hamiltonian method yields A 3 P 0 -1 S 0 = 0.0197 s -1 ). In the process of this calculation, we found out by comparison with the results of mdfgme code [START_REF] Indelicato | Mcdfgme, a multiconfiguration dirac fock and general matrix elements program[END_REF] that non-orthogonality between spinorbitals plays very important role in the evaluation of non-diagonal hyperfine matrix elements (neglecting it may even lead to a sign change). We employed a new code, developed by one of us (P.J.) to evaluate correlated off-diagonal hyperfine matrix elements from GRASP wavefunctions, taking into account the effect of non-orthogonality between spinorbitals.

The accuracy is limited by the electric dipole matrix elements 7s7p 3 P o 1 Q

(1) 1

7s 2 1 S 0 and 7s7p 1 P o 1 Q

(1) 1

7s 2 1 S 0 as well as by the off-diagonal magnetic dipole constants A HF S M1 ( 3 P 1 , 1 P 1 ), A HF S M1 ( 3 P 1 , 3 P 0 ) and A HF S M1 ( 1 P 1 , 3 P 0 ). In case of the matrix element 7s7p 3 P o 1 Q

(1) 1

7s 2 1 S 0 we may assume the 5% relative accuracy for the 3 P 1 -1 S 0 transition rate from experiment [5]. The accuracy of the A( 1 P 1 -1 S 0 ) rate is more difficult to estimate (see ref. [25]). This is the strongest, 'allowed' transition in the radium spectrum, but (as mentioned above) the contribution of its matrix element 7s7p 1 P o 1 Q

(1) 1

7s 2 1 S 0 to the total value of the calculated rate of the hyperfineinduced transition 3 P 0 -1 S 0 is 3.7 times smaller than that of the 7s7p 3 P o 1 Q

(1) 1

7s 2 1 S 0 matrix element, so even a relatively large error bar would be quenched. We may very conservatively take the entire difference between the Babushkin and Coulomb gauge final values from table 1 in reference [25] as the error limit, obtaining a relative accuracy 25% for the A( 1 P 1 -1 S 0 ) transition probability. The accuracy of the diagonal magnetic dipole constants was estimated to be 6% (see [26]). Since an off-diagonal constant depends on both states rather than on one, we, again quite conservatively, doubled the 'diagonal' error limit estimate and assumed 12% as the contributions of the off-diagonal hyperfine constants to the error bar. Eventually, the above procedure yields the lifetime of the metastable 3 P 0 state τ = 34(15)s of 223 88 Ra isotope. The lifetime is based on both the 2-photon and hyperfine-induced channels. The error bar does not include the 2-photon contribution.

The 7s7p 3 P 1 state

Figure 1 presents the transition probability A( 3 P 1 -1 S 0 ) as a function of the size of the multiconfiguration expansion. The transition rates were calculated in Babushkin and Coulomb gauges, which in the non-relativistic formulation correspond to length and velocity form of the transition integral, respectively. Both curves were obtained in the 'core-valence' approximation described above. The resulting Babushkin and Coulomb values are compared with the experimental result obtained by the Argonne group [5] and with three available theoretical values (the value A=4 • 10 6 s -1 obtained by Bruneau [32] did not fit within the vertical scale of figure 1). It is clearly seen in the figure 1 that the core-valence correlation effects are saturated in the framework of the single and restricted double substitutions and five layers of virtual orbitals, as described in section (3) above. The remaining difference between the final Babushkin and Coulomb gauge values may be attributed to the omitted core-core effects. We have made an attempt to estimate the contribution of the core-core correlation to 3 D 2 -1 S 0 and 1 P 1 -1 S 0 transition rates [25], but for other transitions the gauge differences and comparisons with data obtained by other authors, where available, are the only indications of the accuracy of our calculated rates. Although the gauge difference must not be treated as the error bar per se, it is a useful indicator in partially saturated multiconfiguration calculations of transition rates. The values obtained in the Babushkin gauge are weighted toward the outer parts of the electronic wavefunctions, while the Coulomb gauge values weight more inner parts, where the core-core effects arise. Therefore partially saturated expansions often produce Babushkin and Coulomb gauge values converging toward different limits, as in figure 1, with the difference arising from the omitted core-core effects. This is also the reason why Babushkin gauge results are usually treated as more reliable, which seems to be confirmed by the agreement of the Babushkin gauge transition rate A( 3 P 1 -1 S 0 ) with experiment (although agreement this close is very likely accidental; and good agreement for one level is not enough to justify a more general rule). 4.3 The 7s6d 3 D 2 and 7s7p 1 P 1 states

These two levels are distinguished because they were the subject of a separate paper [25]. In the present paper we duly quote the data obtained in core-valence approximation, and the reader is referred to the above mentioned paper for further discussion.

4.4

The remaining 7s7p 3 P 2 , 7s6d 3 D 1 , 7s6d 3 D 3 , and 7s6d 1 D 2 states

The remaining 7s7p 3 P 2 , 7s6d 3 D 1 , 7s6d 3 D 3 , and 7s6d 1 D 2 states, together with the states discussed in sections 4.1, 4.2, and 4.3 above, constitute full set of states arising from the three lowest electronic configurations of radium. Some of these levels should also be considered metastable, not only 7s7p 3 P 0 . The 3 D 3 -3 D 2 transition is the strongest 'direct' decay channel for the 3 D 3 state, but it is in fact very weak (comparable to that of the 3 P 0 -1 S 0 transition). The rates of other possible 'direct' decay channels would be still smaller (e.g. M2 transition 3 D 3 -3 P 1 rate is comparable to that of the 3 D 2 -3 P 0 transition), therefore the multiphoton or hyperfine-induced transitions may also play a role in the 3 D 3 lifetime. Similar considerations may in principle apply to the 3 P 2 , 1 D 2 , and 3 D 2 states (we did not pursue this issue, though). 1 were not used to obtain the lifetimes. With the exception of the 3 P 0 level, the transition rates are given with 4 significant digits and the lifetimes with 3 significant digits, but that not necessarily reflects their accuracy. [28] us [25] to suggest an experimental verification of radium data in Moore's tables [START_REF] Moore | Atomic Energy Levels, NSRDS-NBS No[END_REF]. At that time we were not aware of the papers by Eliav et al [23] and Landau et al [24], where the excitation energies of radium were calculated in the coupled-cluster approximation. More recently these calculations were independently confirmed within the framework of the CI+MBPT theory [28]. In both cases good agreement with experiment had been achieved. There are two principal differences between the methods and approximations used in the abovementioned three papers with respect to the methods and approximations used in the present paper. The transition rates were the primary targets of the present calculations, not the transition energies. Therefore we optimised the electronic wavefunctions separately for each of the nine states of interest. The core-valence correlation effects were fully accounted for. The core-core correlation effects, which are less important in calculations of transition rates and hyperfine structures, were treated in a very crude approximation in the cases of 1 S 0 , 3 D 2 and 1 P 1 states, and were neglected for other states. The calculations of energy level differences require well balanced orbital sets and extensive multiconfiguration expansions. The results are usually in better agreement with experiment if a common set of orbitals is used for both states. If the orbital sets are separately optimised, the transition energy is obtained as a pure difference of the total energies of the two states of interest. They are both several orders of magnitude (five orders in case of radium) larger than the transition energy itself, therefore our calculated transition energy values may be less accurate than the calculated transition rates and hyperfine structures. Any further refinement of the present calculations would require computer resources, which are currently unavailable.

Summary and conclusions

Fig. 1 .

 1 Fig.1. Transition probability A( 3 P1 -1 S0) in Babushkin (upper curve with circles) and Coulomb (lower curve with squares) gauges, as a function of the multiconfiguration expansion (in core-valence approximation -see text for details). The lone dot with error bars (to the right from the end of the Babushkin curve) represents the experimental result from reference[5]. The three stars at the far right represent the theoretical data obtained by other authors and are denoted by publication year in brackets:[1978] -reference[31]; [2000] -reference[START_REF] Dzuba | [END_REF]; [2006] -reference[28].

Table 1

 1 presents calculated transition rates, and the lifetimes are presented in table 2. Coulomb gauge values of electric multipole transitions from the table

Table 1 .

 1 Calculated transition rates between nine lowest levels of radium [s -1 ]. Transition multipolarities are denoted by E1, E2, M1, M2; the Babushkin and Coulomb gauge values by (B) and (C), respectively. HFS means hyperfine-induced. Numbers in brackets represent powers of 10.

	transition	This work Ref. [28] Expt. [5]

Table 2 .

 2 Calculated lifetimes of eight lowest excited states of radium, compared with data from other authors. D3 -3 D2 magnetic dipole transition only; other possible decay channels neglected (c) calculated from transition probabilities quoted after Table VI in Ref.

	state This work Ref. [28]	Ref. [45] Ref. [31] Ref. [32] Expt. [5]
	3 P0	34(15) s a
	3 P1	421 ns	362 ns	505 ns	420 ns	250 ns	422(20) ns
	3 P2	6.46 µs	5.55 µs	5.2 µs
	1 P1	5.56 ns	5.53 ns	5.5 ns
	3 D1	719 µs	654 µs	617 µs
	3 D2	3.95 s	3.3 s	15 s
	3 D3	157 s b
	1 D2	1.37 ms	0.129 ms c 38 ms
	(a) based on E1M1 and hfs-induced decay channels in the 223 88 Ra isotope
	(b) based on 3

Acknowledgments

This work was supported by the Polish Ministry of Science and Higher Education (MNiSW) in the framework of the scientific grant No. 1 P03B 110 30 awarded for the years 2006-2009. Laboratoire Kastler Brossel is Unité Mixte de Recherche du CNRS n • C8552. P.J. acknowledges the support from the Swedish Research Council (VR).

4.025 [3] 3

As discussed in section (4.2) above, the accuracy of our results is difficult to estimate, but it is probably much worse than 3 or 4 digits, particularly for very weak transitions.

The results of the present calculations may be considered as fully converged in core-valence approximation (convergence has indeed been observed for all transitions, similarly to that shown in figure 1), with core-core effects omitted.

With the exception of the 1 P 1 -3 D 2 and 1 D 2 -3 P 1 transitions (we cannot offer any plausible explanation for these discrepancies), our Babushkin gauge values are in reasonably good agreement with the results of Dzuba and Ginges [28]. An interesting feature is the large discrepancy between the results obtained in Babushkin and Coulomb gauges for transitions connecting closely-lying triplet P and triplet D states. The B/C ratios (i.e. the ratios of Babushkin versus Coulomb gauge results) turned out to be much closer to unity for those transitions when the experimental energies in the transition operator were replaced by theoretical ones. This observation, together with inconsistent results of the calculations of energies of excited states, had led