text:
# 1 Introduction

## 1 Introduction

Although there is a considerable literature dealing with Grassmannian spaces (see for example , , , , , , ), the problem of finding the best packings in such spaces seems to have received little attention.

We have made extensive computations on this problem, and have found a number of putatively optimal packings. These computations have led us to conclude that the best definition of distance on Grassmannian space is the “chordal distance” defined in Section 2.

Sections 3, 4, 5, 6 discuss the problems of packing lines in $`R^3`$, planes in $`R^4`$, $`n`$-spaces in $`R^m`$, and lines in $`R^m`$, respectively.

The complete results have been placed in a database on NJAS’s home page. Our search has concentrated on packings of $`N55`$ subspaces of $`G(m,n)`$, for $`m16`$, $`n3`$.

## 2 Grassmannian space

The Grassmannian space $`G(m,n)`$ is the set of all $`n`$-dimensional subspaces of real Euclidean $`m`$-dimensional space $`R^m`$. This is a homogeneous space isomorphic to $`O(m)/(O(n)\times O(mn))`$, and forms a compact Riemannian manifold of dimension $`n(mn)`$.

We first discuss how to define the distance between two $`n`$-planes $`P`$, $`QG(m,n)`$. The principal angles $`\theta _1,\mathrm{},\theta _n[0,\pi /2]`$ between $`P`$ and $`Q`$ are defined by (we follow , p. 584)

$$\mathrm{cos}\theta _i=\underset{uP}{\mathrm{max}}\underset{vQ}{\mathrm{max}}uv=u_iv_i,$$

for $`i=1,\mathrm{},n`$, subject to $`uu=vv=1`$, $`uu_j=0`$, $`vv_j=0`$ $`(1ji1)`$. The vectors $`\{u_i\}`$ and $`\{v_j\}`$ are principal vectors corresponding to the pair $`P`$ and $`Q`$.

Wong  shows that the geodesic distance on $`G(m,n)`$ between $`P`$ and $`Q`$ is

$$d_g(P,Q)=\sqrt{\theta _1^2+\mathrm{}+\theta _n^2}.$$
(2.1)
However, this definition has one drawback: it is not everywhere differentiable. Consider the case $`n=1`$, for example, and hold one line $`P`$ fixed while rotating another line $`Q`$ (both passing through the origin). As the angle $`\varphi `$ between $`P`$ and $`Q`$ increases from 0 to $`\pi `$, the principal angle $`\theta _1`$ increases from 0 to $`\pi /2`$ and then falls to 0, and is non-differentiable at $`\pi /2`$ (see Fig. 1).

Although one might expect this non-differentiability to be a mere technicality, it does in fact cause considerable difficulties for our optimizer, especially in higher dimensions in cases when many distances fall in the neighborhood of singular points of $`d_g`$.

An alternative measure of distance, which we call the chordal distance, is given by

$$d_c(P,Q)=\sqrt{\mathrm{sin}^2\theta _1+\mathrm{}+\mathrm{sin}^2\theta _n}.$$
(2.2)
The reason for the name will be revealed later. This approximates the geodesic distance when the planes are close, has the property that its square is differentiable everywhere, and, as we shall attempt to demonstrate, has a number of other desirable features.

A third definition has been used by Asimov  and Golub and Van Loan , p. 584, namely

$$d_m(P,Q)=\theta _n=\underset{i=1,\mathrm{},n}{\mathrm{max}}\theta _i.$$

This shares the vices of the geodesic distance.

Of course for $`n=1`$ all three definitions are equivalent, in the sense that they lead to the same optimal packings.

We can now state the packing problem: given $`N,n,m`$, find a set of $`n`$-planes $`P_1,\mathrm{},P_NG(m,n)`$ so that $`\underset{ij}{\mathrm{min}}d(P_i,P_j)`$ is as large as possible, where $`d`$ is either geodesic or chordal distance. Since $`G(m,n)`$ is compact, the problem is well-defined. Because $`G(m,n)`$ and $`G(m,mn)`$ are essentially the same space, we will usually assume that $`nm/2`$.

We also need some further terminology. A generator matrix for an $`n`$-plane $`PG(m,n)`$ is an $`n\times m`$ matrix whose rows span $`P`$. The orthogonal group $`O(m)`$ acts on $`G(m,n)`$ by right multiplication of generator matrices. The automorphism group of a subset $`\{P_1,\mathrm{},P_N\}G(m,n)`$ is the subset of $`O(m)`$ which fixes or permutes these planes.

By applying a suitable element of $`O(m)`$ and choosing appropriate basis vectors for the planes, any given pair of $`n`$-planes $`P,Q`$ with $`nm/2`$ can be assumed to have generator matrices

$$\left[\begin{array}{ccccccccccc}1& 0& \mathrm{}& 0& 0& 0& \mathrm{}& 0& 0& \mathrm{}& 0\\ 0& 1& \mathrm{}& 0& 0& 0& \mathrm{}& 0& 0& \mathrm{}& 0\\ & & \mathrm{}& & & & \mathrm{}& & & \mathrm{}& \\ 0& 0& \mathrm{}& 1& 0& 0& \mathrm{}& 0& 0& \mathrm{}& 0\end{array}\right]$$
(2.3)
and

$$\left[\begin{array}{ccccccccccc}\mathrm{cos}\theta _1& 0& \mathrm{}& 0& \mathrm{sin}\theta _1& 0& \mathrm{}& 0& 0& \mathrm{}& 0\\ 0& \mathrm{cos}\theta _2& \mathrm{}& 0& 0& \mathrm{sin}\theta _2& \mathrm{}& 0& 0& \mathrm{}& 0\\ & & \mathrm{}& & & & \mathrm{}& & & \mathrm{}& \\ 0& 0& \mathrm{}& \mathrm{cos}\theta _n& 0& 0& \mathrm{}& \mathrm{sin}\theta _n& 0& \mathrm{}& 0\end{array}\right]$$
(2.4)
respectively, where $`\theta _1,\mathrm{},\theta _n`$ are the principal angles between them (, Theorem 2).

## 3 Packing lines through the origin in 3-space

Our initial work on this problem was prompted by a question raised in 1992 by Julian Rosenman, an oncologist at the Univ. of North Carolina School of Medicine and Computer Science, in connection with the treatment of tumors using high energy laser beams . He asked for the best way to separate $`N`$ lines through a given point in $`R^3`$, or in other words for the best packings in $`G(3,1)`$.

Together with W. D. Smith, we had been carrying out an extensive search for the best packings of a given number of points on $`S^2`$, i.e. spherical codes , , and we therefore modified our programs to search instead for packings of lines. We omit the details of this search, since we later recalculated these results using the more general methods described in Section 4.

The results are summarized in Table 1, which gives, for each value of $`N`$ in the range 2 to 55, the minimal angle of the best packing we have found of $`N`$ lines, and for comparison the minimal angle of the best packing known of $`2N`$ points on $`S^2`$ (taken from ). We see that requiring a packing of $`2N`$ points on $`S^2`$ to be antipodal is a definite handicap: only in the cases of 6 and 12 points do the antipodal and unrestricted packings coincide. Decimals in the tables have been rounded to four decimal places.

The last two columns of the table specify the largest automorphism group111That is, the subgroup of $`O(3)`$ that fixes or permutes the $`2N`$ points. we have found of any such best antipodal packing of $`2N`$ points. The fourth column gives the order of the group and the fifth and sixth columns its name in the orbifold notation (cf. ) and as the double cover of a rotation group. The symbol $`\pm 𝒢`$ indicates that the group consists of the matrices $`\pm M`$ for $`M𝒢`$, where $`𝒢`$ is a cyclic $`(𝒞)`$, dihedral $`(𝒟)`$, tetrahedral $`(𝒯)`$, octahedral $`(𝒪)`$ or icosahedral $`()`$ group. In each case the subscript gives the order of the rotation group.

In some cases the best packings can be obtained by taking the diameters of a known polyhedron, and if so this is indicated in the final column of the table. Other entries in this column refer to the brief descriptions given below.

The entries for $`N6`$ were shown to be optimal by Fejes Tóth in 1965  (see also Rosenfeld ), and the 7-line arrangement will be proved optimal in Section 5. The solutions for $`N8`$ are the best found with over 15000 random starts with our optimizer. There is no guarantee that these are optimal, but experience with similar problems suggests that they will be hard to beat and in any case will be not far from optimal.

For $`N=1,2,3,6`$ the solutions are known to be unique, for $`N=4`$ there are precisely two solutions (, , ), and for $`N=5,7,8`$ the solutions appear to be unique. For larger values of $`N`$, however, the solutions are often not unique. For $`N=9`$ lines there are two different solutions, and in the range $`N30`$ the solutions for 10, 22, 25, 27, 29 lines (and possibly others) contain lines that “rattle”, that is, lines which can be moved freely over a small range of angles without affecting the minimal angle.

### Notes on Table 1

* Five of the six diameters of a regular icosahedron.
* The putatively optimal arrangement forms an unpleasant-looking configuration of 16 antipodal points on $`S^2`$, with no further symmetry. The convex hull is shown in Fig. 2. In contrast, the putatively best packing of 16 points (Fig. 3) has a group of order 16.

#### $`N=9`$.

There appear to be two inequivalent solutions. The nicest has symmetry group $`3\times `$, of order 12 (or $`[2^+,6^+]`$ in Coxeter’s notation ) — see Fig. 4. The points lie in six horizontal layers of equilateral triangles. The nine points in the Northern hemisphere can be located by drawing seven equilateral spherical triangles of edge length 47.98213264, arranged as in Fig. 5. The North pole is located at the midpoint of the central triangle. The second solution, shown in Fig. 6, has group $`2`$, of order 4.

#### $`N=10`$.

This is a “hexakis bi-antiprism”, since it consists of a bi-(hexagonal antiprism) together with an axial line that can “rattle”, giving infinitely many solutions. The axis appears horizontally in Fig. 7.

#### $`N=15`$.

Combinatorially, this is an “axially depleted pentakis dodecahedron”, obtained from a pentakis dodecahedron (the solution for $`N=16`$) by omitting two opposite vertices (see Fig. 8). However, the angular separation is slightly greater than for $`N=16`$.

#### $`N=55`$.

The 110 antipodal points can be taken to be the union of the vertex sets of a dodecahedron (20), an icosidodecahedron (30) and a truncated icosahedron (60) — see Fig. 9. The 15 lines through the icosidodecahedral points “rattle”. This arrangement of points can be obtained from a “geodesic dome” or 12-fold reticulated icosahedron  by omitting its twelve pentagonal vertices, and so may be described as consisting of the vertices of a “depleted 12-fold reticulated icosahedron”. Incidentally, our best solution for 61 lines has a symmetry group of order 6, and is not found by placing points at the centers of the pentagonal faces of Fig. 9.

Packings of lines in higher dimensional spaces will be discussed in Section 6.

## 4 Packing two-dimensional planes in $`R^4`$

In  Asimov proposed a technique called the “Grand Tour” for displaying multi-dimensional data on a two-dimensional computer screen. His idea was to choose a finite sequence or “tour” of two-dimensional planes that are in some sense “representative” of all planes, and to project the data onto each plane in turn, in the hope that the viewer will be able to notice any pattern or structure that is present. This technique has been implemented in the XGobi program , .

In 1993 Dianne Cook asked us if we could modify our algorithm for finding spherical codes, in order to search for packings in $`G(m,2)`$ for $`m4`$. Furthermore, for the Grand Tour application, there should be a Hamiltonian circuit through the planes, in which edges indicate neighboring planes. We accomplished this in the following way.

Following the methods we had used to find packings in $`S^d`$ and experimental designs in various spaces (see , Sec. 4.11; , ), corresponding to a set of planes $`S=\{P_1,\mathrm{},P_N\}G(m,n)`$ we define a potential

$$\mathrm{\Phi }_c(S)=\underset{i<j}{}\frac{1}{d_c(P_i,P_j)A},$$

where $`A`$ is a suitably chosen constant. There is a similar definition for $`\mathrm{\Phi }_g`$ involving $`d_g`$. Initially $`A`$ is set to 0 and $`S`$ is a randomly chosen set of planes. We invoke the Hooke-Jeeves “pattern search” optimization algorithm to modify $`S`$, attempting to minimize $`\mathrm{\Phi }_c(S)`$. After a fixed time (we used 100 steps of the optimizer), $`A`$ is advanced halfway to the current minimal distance between the $`P_i`$, and the process is repeated, terminating when no further improvement is obtained to the accuracy of the machine. The whole process is repeated with several thousand random starts — and also with initial configurations taken from other sources, such as subsets or supersets of other arrangements — and the best final configuration is recorded.

The partial derivatives of $`\mathrm{\Phi }_c(S)`$ were found analytically, but for $`\mathrm{\Phi }_g(S)`$ we calculated them by numerical differentiation.

We began by computing a table of packings of $`N`$ planes in $`G(4,2)`$, for $`N50`$, with respect to both $`d_c`$ and $`d_g`$. The results for packings of $`N24`$ planes are summarized in Table 2. Columns 2 and 3 refer to the best packings found for the chordal distance, and give both

$$d_c^2=\underset{ij}{\mathrm{min}}d_c^2(P_i,P_j)\mathrm{and}d_g^2=\underset{ij}{\mathrm{min}}d_g^2(P_i,P_j)$$

for such packings. Columns 4 and 5 refer similarly to the best packings found with respect to geodesic distance.

Understanding these results was hindered by the fact that planes are, as the name suggests, plain, with no distinguishing features, and when produced by the computer appear as random generator matrices referred to a random coordinate frame. However, we found that the set of principal vectors in the planes could often be used to find a coordinate system that would reveal some of the structure of the planes.

For two planes the best packing is the same for both metrics, and consists of two mutually orthogonal planes with principal angles $`\pi /2`$, $`\pi /2`$, so $`d_c^2=2`$, $`d_g^2=\pi ^2/2`$.

For three planes the two answers are different. The best packing for chordal distance (Fig. 10a) has principal angles $`\pi /3`$, $`\pi /3`$ between each pair of planes, so $`d_c^2=3/2`$, $`d_g^2=2\pi ^2/9`$, whereas the best packing for geodesic distance (Fig. 10b) has angles $`\pi /6`$, $`\pi /2`$, so $`d_c^2=5/4`$, $`d_g^2=5\pi ^2/18`$.

Postponing discussion of $`N=4`$ and 5 for the moment, let us consider the case of six planes, where we discovered that the answer for both metrics formed a regular simplex222Although not an isoclinic configuration, cf. Wong . There seems to be only a slight overlap connection between that problem and ours.: the principal angles between every pair of planes were arcsin $`1/\sqrt{5}`$ and $`\pi /2`$. This set of planes is conveniently described using simplicial coordinates. Let $`A`$, $`B`$, $`C`$, $`D`$, $`E`$ be the vectors from the center of a regular simplex in $`R^4`$ to its vertices, with $`A+B+C+D+E=0`$, and write $`[a,b,c,d,e]`$ for $`aA+bB+cC+dD+eE`$. Then one of the six planes is spanned by $`[1,\tau ,1,0,0]`$ and its cyclic shifts, where $`\tau =(1+\sqrt{5})/2`$, and the other five planes are obtained from it by taking even permutations of these coordinates. (We later found an equally good packing with respect to chordal distance that was not a simplex with respect to geodesic distance. This is described below.)

The six planes intersect the surface of the unit ball of $`R^4`$ in a remarkable link, the “Hexalink”. It can be shown that apart from the $`n`$-component unlinks and Hopf links, which exist for all $`n`$, the Hexalink is the only link of circular rings in which there are orientation-preserving symmetries taking any two links to any other two. (Details will be given elsewhere.)

The discovery of this 6-vertex simplex was initially somewhat of a surprise, since $`G(4,2)`$ is only a 4-dimensional manifold. It did suggest that the answer should somehow be related to the six equi-angular diameters of the icosahedron, and led to the following reformulation of the problem.

We remind the reader that any element $`\alpha `$ of $`SO(4)`$ may be represented as

$$\alpha :x\overline{\mathrm{}}xr,$$

where $`x=x_0+x_1i+x_2j+x_3k`$ represents a point on $`S^3`$ and $`\mathrm{}`$, $`r`$ are unit quaternions . The pair $`\mathrm{}`$, $`r`$ represent the same $`\alpha `$. The correspondence between $`\alpha `$ and $`\pm (\mathrm{},r)`$ is one-to-one.

Given a plane $`PG(4,2)`$, let $`\alpha `$ be the element of $`SO(4)`$ that fixes $`P`$ and negates the points of the orthogonal plane $`P^{}`$. Then $`\alpha ^2=1`$, and for this $`\alpha `$, it is easy to see that $`\mathrm{}=\mathrm{}_1i+\mathrm{}_2j+\mathrm{}_3k`$ and $`r=r_1i+r_2j+r_3k`$ are purely imaginary unit quaternions. This establishes the following result (which can be found for example in Leichtweiss ).

###### Theorem 1

A plane $`PG(4,2)`$ is represented by a pair $`(\mathrm{},r)S^2\times S^2`$, with $`(\mathrm{},r)`$ representing the same $`P`$.

There are simple formulae relating $`P`$ and $`(\mathrm{},r)`$, pointed out to us by Simon Kochen. Given $`(\mathrm{},r)`$, if $`\mathrm{}r`$ then $`P`$ is spanned by the vectors corresponding to the quaternions $`u=1\mathrm{}r`$ and $`v=\mathrm{}+r`$. The special case when $`\mathrm{}=r`$ is even simpler: take $`u=1`$, $`v=\mathrm{}`$. If $`\mathrm{}=r`$, take $`u`$ and $`v`$ to be purely imaginary unit quaternions such that $`\mathrm{},u,v`$ correspond to a coordinate frame. Conversely, if $`P`$ is spanned by two orthogonal unit vectors represented by quaternions $`u`$, $`v`$, then $`\pm (\mathrm{},r)=\pm (u\overline{v}v\overline{u},\overline{v}u\overline{u}v)`$.

Given two planes $`P,QG(4,2)`$, represented by $`\pm (\mathrm{},r)`$, $`\pm (\mathrm{}^{},r^{})`$, respectively, the principal angles $`\theta _1`$, $`\theta _2`$ between them may be found as follows. Let $`\varphi `$ (resp. $`\psi `$) be the angle between $`\mathrm{}`$ and $`\mathrm{}^{}`$ (resp. $`r`$ and $`r^{}`$), with $`0\varphi ,\psi \pi `$. If $`\varphi +\psi >\pi `$, replace $`\varphi `$ by $`\pi \varphi `$ and $`\psi `$ by $`\pi \psi `$, so that $`0\varphi +\psi \pi `$, with $`\varphi \psi `$ (say). Then

$$\theta _1,\theta _2=\frac{\psi \pm \varphi }{2},$$
(4.5)
$`d_g^2(P,Q)`$ $`=`$ $`{\displaystyle \frac{\psi ^2+\varphi ^2}{2}},`$ (4.6)
$`d_c^2(P,Q)`$ $`=`$ $`1\mathrm{cos}\psi \mathrm{cos}\varphi .`$ (4.7)
We omit the elementary proof.

A set $`S=\{P_1,\mathrm{},P_N\}G(4,2)`$ is thus represented by a “binocular code” consisting of a set of pairs $`\pm (\mathrm{}_i,r_i)S^2\times S^2`$. We call the list of $`2N`$ points $`\pm \mathrm{}_i`$ (they need not be distinct) the “left code” corresponding to $`S`$, and the points $`\pm r_i`$ the “right code”. Conversely, given two multisets $`LS^2`$, $`RS^2`$, each of size $`2N`$ and closed under negation, and a bijection or “matching” $`f`$ between them that satisfies $`f(\mathrm{})=f(\mathrm{})`$, $`\mathrm{}L`$, we obtain a set of $`N`$ planes in $`G(4,2)`$.

The binocular codes for the $`d_c`$-optimal packings of $`N=2,\mathrm{},6`$ planes are shown in Figs. 11, 12. Except for $`N=3`$, the left and right codes are identical. Matching points from the left and right codes are labeled with the same symbol. For $`N5`$ the points lie in the equatorial plane, and for $`N4`$ there are repeated points. The points lie on regular figures, except for $`N=4`$ where the points are $`\pm (1,0,0)`$, $`\pm (\frac{1}{\sqrt{3}},\pm \sqrt{\frac{2}{3}},0)`$.

For $`N=6`$ the left and right codes consist of the 12 vertices of an icosahedron (Fig. 12). Let these be the points

$$\lambda (0,\pm 1,\pm \tau ),\lambda (\pm \tau ,0,\pm 1),\lambda (\pm 1,\pm \tau ,0),$$
(4.8)
where $`\lambda =1/\sqrt{\tau +2}`$. The matching is obtained by mapping each point to its algebraic conjugate (i.e. replacing $`\sqrt{5}`$ by $`\sqrt{5}`$), and rescaling so the points again lie on a unit sphere. As already mentioned, the principal angles between each pair of these planes are $`\mathrm{arcsin}1/\sqrt{5}`$ and $`\pi /2`$, so $`d_c^2=6/5`$, $`d_g^2=2.6824`$.

There is a second set of six planes with $`d_c^2=6/5`$, but with $`d_g^2`$ only equal to 2.0030. Here the left and right codes form the vertices of what we shall call the “anti-icosahedron”, consisting of the points $`\lambda (0,\pm 1,\pm \tau )`$, $`\lambda (\pm \tau ,0,\pm 1)`$, $`\lambda (\pm \tau ,\pm 1,0)`$ (compare (4.8)), and shown in Fig. 13. Topologically this is a “parallel bi-slit cuboctahedron”, obtained by dividing two opposite square faces of a cuboctahedron into two triangles by parallel lines. The matching sends $`\lambda (0,\pm 1,\tau )`$ to $`\lambda (0,\pm \tau ,1)`$, $`\lambda (\tau ,0,\pm 1)`$ to $`\lambda (1,0,\pm \tau )`$ and $`\lambda (\pm \tau ,1,0)`$ to $`\lambda (\pm 1,\tau ,0)`$.

Four larger putatively $`d_c`$-optimal packings that the computer found are also worth mentioning, those of 10, 18, 48 and 50 planes.

For the 10-plane packing the left and right codes consist of the vertices of a decagonal prism with coordinates

$$\pm (\sqrt{\frac{2}{3}}\mathrm{cos}r\theta ,\sqrt{\frac{2}{3}}\mathrm{sin}r\theta ,\sqrt{\frac{1}{3}}),r=0,\mathrm{},9,$$

where $`\theta =\pi /5`$. A typical point $`(\sqrt{2/3}\mathrm{cos}r\theta ,\sqrt{2/3}\mathrm{sin}r\theta ,\sqrt{1/3})`$ in the left code is matched with $`(1)^r`$ $`(\sqrt{2/3}\mathrm{cos}\mathrm{\hspace{0.17em}3}r\theta ,\sqrt{2/3}\mathrm{sin}3r\theta ,\sqrt{1/3})`$ in the right code. Then $`d_c^2=10/9`$, $`d_g^2=1.5725`$. The most interesting property of this configuration is that, although three different sets of canonical angles occur, the chordal distance between every pair of planes is the same — this is a regular simplex!

We find it surprising that this is superior to the packing of ten planes obtained by matching the vertices of a dodecahedron to their algebraic conjugates, as shown in Figure 14.

#### 18 planes.

Let $`𝒪`$ denote the set of vertices of a regular octahedron. The binocular code for the 18-plane packing is $`𝒪\times 𝒪:=\{(\mathrm{},r):\mathrm{},r𝒪\}`$. Projectively, the 18 two-spaces are the edges of a desmic triple333Two tetrahedra are called desmic if they are in perspective from four distinct points — these form a third tetrahedron desmic with either. The 18 edges of such a desmic triple are also the edges of a unique other desmic triple. of tetrahedra , , . We were pleased to see this configuration appear, since it was already familiar in the context of quantum logic . This set of planes can be constructed directly by taking the planes spanned by all pairs of mutually perpendicular minimal vectors of the $`D_4`$ lattice (cf. ). In this form the planes have generator matrices such as

$$\left[\begin{array}{cccc}+& +& 0& 0\\ 0& 0& +& \end{array}\right](12\mathrm{planes}),\left[\begin{array}{cccc}+& +& 0& 0\\ +& & 0& 0\end{array}\right](6\mathrm{planes}).$$

Three different sets of principal angles occur, $`0`$, $`\pi /2`$; $`\pi /4`$, $`\pi /4`$; and $`\pi /2`$, $`\pi /2`$; so that $`d_c^2=1`$, $`d_g^2=1.2337`$. The automorphism group of this packing has structure $`[3,4,3].2`$ and order 2304. These 18 planes can also be described as the set of planes that meet the 24-cell in one of its equatorial squares.

#### 48 and 50 planes.

Let $`𝒞`$ denote the set of vertices of a cube. Then the binocular codes for the 48- and 50-plane packings are respectively $`𝒪\times 𝒞𝒞\times 𝒪`$ and $`𝒪\times 𝒪𝒞\times 𝒞`$. For both arrangements we have $`d_c^2=2/3`$, $`d_g^2=0.7576`$. Both are believed to be optimal. (We mention the 48-plane packing because of its greater symmetry.)

#### The matching problem.

In order to maximize the minimal chordal distance between the planes, the matching between the left and right codes should (from (4.7)) be chosen so as to minimize the maximal value of $`\mathrm{cos}\psi \mathrm{cos}\varphi `$. Stated informally, the matching should be such that if two points are close together on one of the spheres then the points to which they are matched should be far apart on the other sphere. In the case of the icosahedron, for example (see Fig. 12), the matching sends adjacent vertices to non-adjacent vertices. There is a unique way to do this.

At this point we were tempted to see if any new record packings in $`G(4,2)`$ could be obtained by taking the $`2N`$ antipodal points corresponding to a good $`N`$-line packing in $`G(3,1)`$, and matching them with themselves in an optimal way. However, it was not easy to see how to solve the matching problem. Fortunately David Applegate (personal communication) found that it could be reformulated as an integer programming problem, as follows.

Given an antipodal set $`S=\{P_1,\mathrm{},P_{2N}\}S^2`$, we wish to find a permutation $`f`$ of $`S`$ with the property that $`f(P_i)=f(P_i)`$, for all $`i`$, and such that the minimal value of $`1P_iP_jf(P_i)f(P_j)`$ $`(ij)`$ is maximized. If we do the maximization by binary search, we may define $`\pi (P_i,P_j)=1`$ if $`f(P_i)=P_j`$, or 0 otherwise, with the constraints

$`\pi (P_i,P_j)`$ $`=`$ $`\pi (P_i,P_j),`$
$`{\displaystyle \underset{j=1}{\overset{2N}{}}}\pi (P_i,P_j)`$ $`=`$ $`1,1i2N,`$
$`{\displaystyle \underset{i=1}{\overset{2N}{}}}\pi (P_i,P_j)`$ $`=`$ $`1,1j2N,`$
and

$$\pi (P_i,P_j)+\pi (P_k,P_{\mathrm{}})<1\text{for all}1i,j,k,\mathrm{}2N\text{such that }1P_iP_kP_jP_{\mathrm{}}<M.$$

There is a feasible solution if and only if $`d_c^2M`$.

Applegate kindly implemented this procedure, and used it to solve the matching problem for our best packings in $`G(3,1)`$ and for various polyhedra. Unfortunately no new records have yet been obtained by this method.

#### Best packings for the geodesic metric.

So far we have mostly discussed packings that attempt to maximize the minimal chordal distance between planes. We also computed packings for the geodesic distance, and we shall now describe some of them. In general, however, these are much less symmetric than the chordal-distance packings, especially for more than 16 planes.

For $`N=4,5`$ the solutions are subsets of the 6-plane packing. For $`N=7`$ the left code is a heptagonal anti-prism and the right code is an equatorial 14-gon.

For $`N=12`$ the binocular code is the set $`\{\pm (\mathrm{},r)\}`$, where $`\mathrm{}`$ (resp. $`r`$) runs through the vertices of a regular tetrahedron (resp. equatorial equilateral triangle).

For $`N=16`$ the left and right codes are the union of two similarly oriented square prisms. For larger $`N`$ the geodesic packings do not seem so interesting. For $`N=18`$, for example, the best packing has no non-trivial symmetries.

#### Hamiltonian paths.

One further question remains to be discussed, that of arranging the planes in a circuit in such a way that adjacent planes are close together, for the “Grand Tour” application. This turned out to be a much easier problem than finding the packings. We handled it in two different ways. For configurations such as the 48-plane arrangement, where there was an obvious notion of adjacency — in this case, defining two planes to be adjacent if the principal angles are $`\pi /4`$, $`\pi /4`$ — we represent the packing by a graph with nodes representing the planes, and look for a Hamiltonian cycle. In less regular cases, we convert the packing into a traveling salesman problem using chordal distance to define the distance between nodes, and look for a minimal length circuit. In both cases we were able to make use of the travelling salesman programs of Applegate et al. , which can handle 100-node graphs without difficulty. Some of the packings in the database mentioned in Sect. 1 (those with suffix .ham) have been arranged in cycles in this way.

The above reformulation in terms of binocular codes applies only to $`G(4,2)`$ (we realized from the beginning that this case would be special, since $`G(4,2)`$ is the only Grassmannian space where the Riemannian metric is not unique ). In the next section we describe a second reformulation that applies to the general case.

## 5 Packing $`n`$-planes in $`R^m`$

Three observations contributed to the second reformulation.

(i) We noticed (see Table 2) that for several examples of $`N`$-plane packings in $`G(4,2)`$ the largest value of $`d_c^2`$ that we could attain was $`N/(N1)`$, and that in every case this was an upper bound. Further experimentation with other packings in $`G(m,n)`$ for $`m8`$ led us to guess an upper bound of

$$d_c^2\frac{n(mn)}{m}\frac{N}{N1},$$
(5.9)
which again we could achieve for some small values of $`N`$. The form of (5.9) was suggestive of the Rankin bound for spherical codes  or the Plotkin bound for binary codes .

(ii) Investigation of the 18-plane chordal-distance packing revealed that, with respect to chordal distance, this has the structure of a regular orthoplex (or generalized octahedron) with 18 vertices. Combining this with the fact that the 10-plane packing formed a regular simplex, we had strong evidence that $`G(4,2)`$ should have an isometric embedding (with respect to $`d_c`$) into $`R^9`$.

(iii) A computer program was therefore written to determine the lowest dimensions into which our library of packings in $`G(m,n)`$ could be isometrically embedded. More precisely, for a given set of $`N`$ planes in $`G(m,n)`$, we searched for the smallest dimension $`D`$ such that there are $`N`$ points in $`R^D`$ whose Euclidean distances coincide with the chordal distances between our planes.

The results were a surprise: it appeared that $`G(m,n)`$ with chordal distance could be isometrically embedded into $`R^D`$, for $`D=\left(\genfrac{}{}{0pt}{}{m+1}{2}\right)1`$, independent of $`n`$. Furthermore the points representing elements of $`G(m,n)`$ were observed to lie on a sphere of radius $`\sqrt{n(mn)/2m}`$ in $`R^D`$.

Aided by discussions with Colin Mallows, we soon found an explanation: just associate to each $`PG(m,n)`$ the orthogonal projection map from $`R^m`$ to $`P`$. If $`A`$ is a generator matrix for $`P`$ whose rows are orthogonal unit vectors, then the projection is represented by the matrix

$$𝒫=A^{tr}A.$$
(5.10)
$`𝒫`$ is an $`m\times m`$ symmetric idempotent matrix, which is independent of the particular orthonormal generator matrix used to define it. Changing to a different coordinate frame in $`R^m`$ has the effect of conjugating $`𝒫`$ by an element of $`O(m)`$. With the help of (2.3), we see that

$$\mathrm{trace}𝒫=n.$$
(5.11)
Thus $`𝒫`$ lies in a space of dimension $`\left(\genfrac{}{}{0pt}{}{m+1}{2}\right)1`$.

Let $``$ denote the $`L_2`$-norm of a matrix: if $`M=(M_{ij})`$, $`1i,jm`$,

$$M=\sqrt{\underset{i=1}{\overset{m}{}}\underset{j=1}{\overset{m}{}}M_{ij}^2}=\sqrt{\mathrm{trace}M^{tr}M}.$$

For $`P,QG(m,n)`$, with orthonormal generator matrices $`A`$, $`B`$, and principal angles $`\theta _1,\mathrm{},\theta _n`$, an elementary calculation using (2.3), (2.4) shows that

$`d_c^2(P,Q)`$ $`=`$ $`n(\mathrm{cos}^2\theta _1+\mathrm{}+\mathrm{cos}^2\theta _n)`$ (5.12)
$`=`$ $`n\mathrm{trace}A^{tr}AB^{tr}B`$
$`=`$ $`{\displaystyle \frac{1}{2}}𝒫𝒬^2,`$
where $`𝒫`$, $`𝒬`$ are the corresponding projection matrices.

Note that if we define the “de-traced” matrix $`\overline{𝒫}=𝒫\frac{n}{m}I_m`$, then $`\mathrm{trace}\overline{𝒫}=0`$, and $`\overline{𝒫}^2=\frac{n(mn)}{n}`$. We have thus established the following theorem.

###### Theorem 2

The representation of $`n`$-planes $`PG(m,n)`$ by their projection matrices $`\overline{𝒫}`$ gives an isometric embedding of $`G(m,n)`$ into a sphere of radius $`\sqrt{n(mn)/n}`$ in $`R^D`$, $`D=\left(\genfrac{}{}{0pt}{}{m+1}{2}\right)1`$, with $`d_c(P,Q)=\frac{1}{\sqrt{2}}𝒫𝒬`$.

Thus chordal distance between planes is $`1/\sqrt{2}`$ times the straight-line distance between the projection matrices (which explains our name for this metric). The geodesic distance between the planes is $`1/\sqrt{2}`$ times the geodesic distance between the projection matrices, measured along the sphere in $`R^D`$.

Figure 5 attempts to display the embeddings of $`G(m,0)`$, $`G(m,1),\mathrm{},G(m,m)`$ in $`R^{D+1}`$. Since $`𝒫\frac{1}{2}I_m^2=m/4`$, all the planes lie on the large sphere, centered at $`\frac{1}{2}I_m`$, of radius $`\sqrt{m}/2`$. The members of $`G(m,n)`$ lie on the intersection of the large sphere with the plane $`\mathrm{trace}(𝒫)=n`$, which intersection is itself a sphere in $`R^D`$ of radius $`\sqrt{n(mn)/m}`$ centered at $`\frac{n}{m}I_m`$. A plane $`P`$ and its orthogonal complement $`P^{}`$ are represented by antipodal points on the large sphere.

In contrast to this result, we briefly remark, without giving details, that there is no way to embed $`G(m,n)`$ into Euclidean space of any dimension so that the geodesic distance $`d_g`$ on $`G(m,n)`$ is represented by Euclidean distance in that space. Of course the Plücker embedding, in which members of $`G(m,n)`$ are represented by points in projective space of dimension $`\left(\genfrac{}{}{0pt}{}{m}{n}\right)1`$, also does not give a way to realize either $`d_c`$ or $`d_g`$ as Euclidean distance. (Nor does the Nash embedding theorem .) Note also that the dimension of the Plücker embedding is in general much larger than the dimension of our embedding.

Since we have embedded $`G(m,n)`$ into a sphere of radius $`\sqrt{n(mn)/m}`$ in $`R^D`$, we can apply the Rankin bounds for spherical codes , and deduce:

#### Corollary.

(i) The simplex bound: for a packing of $`N`$ planes in $`G(m,n)`$,

$$d_c^2\frac{n(mn)}{m}\frac{N}{N1}.$$
(5.13)
Equality requires $`ND+1=\left(\genfrac{}{}{0pt}{}{m+1}{2}\right)`$, and occurs if and only if the $`N`$ points in $`R^D`$ corresponding to the planes form a regular ‘equatorial’ simplex.

(ii) The orthoplex bound: for $`N>\left(\genfrac{}{}{0pt}{}{m+1}{2}\right)`$,

$$d_c^2\frac{n(mn)}{m}.$$
(5.14)
Equality requires $`N2D=(m1)(m+2)`$, and occurs if the $`N`$ points form a subset of the $`2D`$ vertices of a regular orthoplex. If $`N=2D`$ this condition is also necessary.

Lemmons and Seidel (, Theorem 3.6) give a bound for equi-isoclinic packings in $`G(m,n)`$ which agrees with (5.13); of course our bound is more general. The case $`n=1`$ of (5.13) is given in Theorem 3 of .

We were happy to obtain confirmation of (5.9). Since $`d_c^2`$ can never exceed $`n`$, we can also write

$$d_c^2\mathrm{min}\{n,\frac{n(mn)}{m}\frac{N}{N1}\}.$$
(5.15)
The Corollary allows us to establish the optimality of many of our packings. In the range $`m16`$, $`n3`$, $`N55`$, there are over 750 cases which appear to meet (5.15) or (5.14), in the sense that the ratio of $`d_c^2`$ to the bound is greater than .9999999 . For $`n=2`$, for instance, the following are the packings of $`N`$ planes in $`G(m,2)`$ for $`m10`$ that appear to achieve either (5.15) (listed before the semicolon) or (5.14) (after):

$$\begin{array}{cc}\hfill m=4:& 28,10;1118\hfill \\ \hfill m=5:& 411;16,17\hfill \\ \hfill m=6:& 314;22,23\hfill \\ \hfill m=7:& 618;29\hfill \\ \hfill m=8:& 421,28;3744\hfill \\ \hfill m=9:& 625;\hfill \\ \hfill m=10:& 528;\hfill \end{array}$$

However, this is not as meaningful as it at first seems. Consider the packings of $`N`$ planes in $`G(10,2)`$. The ratios of $`d_c^2`$ to the simplex bound for $`N=21,\mathrm{},33`$, rounded to 14 decimal places, are as follows:

$$\begin{array}{cc}21\hfill & 0.99999999999999\hfill \\ 22\hfill & 1.00000000000000\hfill \\ 23\hfill & 0.99999999999999\hfill \\ 24\hfill & 0.99999999999999\hfill \\ 25\hfill & 1.00000000000000\hfill \\ 26\hfill & 0.99999999999999\hfill \\ 27\hfill & 1.00000000000000\hfill \\ 28\hfill & 0.99999999921354\hfill \\ 29\hfill & 0.99999470290071\hfill \\ 30\hfill & 0.99998661084736\hfill \\ 31\hfill & 0.99961159256647\hfill \\ 32\hfill & 0.99909699728979\hfill \\ 33\hfill & 0.99854979246655\hfill \end{array}$$

From this it seems very likely that there exists a 27-plane packing meeting the simplex bound, but certainly further investigation is needed to determine if there is a 28-plane simplex packing in the neighborhood of the computer’s approximate solution. We have carried out such analyses in all the putative cases of equality for packing lines (see Sect. 6), for planes in $`G(4,2)`$ (as already discussed) and for 70 planes in $`G(8,4)`$ (see below). The other cases remain to be investigated.

In principle there is no difficulty in settling these questions. All the chordal distances between the planes are specified. So we could simply set up a (large) system of quartic equations with integer coefficients for the entries in the generator matrices, and ask if there is a real solution.

The chief difficulty in attempting to understand the computer-generated packings is that there are usually very large numbers of equally good solutions. The case of 18 planes in $`G(4,2)`$ is typical in this regard. There is one exceptionally pleasing solution (described in Section 4), in which only three different pairs of principal angles occur. However, it seems that there is a roughly 29-dimensional manifold of solutions, in which a typical solution, although still having $`d_c^2=1`$, has an apparently random set of principal angles. (Our investigation of this question is not yet completed.)

The case of subspaces of dimension $`n`$ in $`R^{2n}`$ or $`R^{2n+1}`$ is especially interesting. The largest possible arrangements of planes that could achieve the two bounds are:

$$\begin{array}{cccccc}m& n& \multicolumn{2}{c}{N(\mathrm{simplex})}& \multicolumn{2}{c}{N(\mathrm{orthoplex})}\\ \hfill 2& \hfill 1& \hfill 3& \hfill & \hfill 4& \hfill \\ \hfill 3& \hfill 1& \hfill 6& \hfill & \hfill 10& \\ \hfill 4& \hfill 2& \hfill 10& \hfill & \hfill 18& \hfill \\ \hfill 5& \hfill 2& \hfill 15& & \hfill 28& \\ \hfill 6& \hfill 3& \hfill 21& & \hfill 40& \\ \hfill 7& \hfill 3& \hfill 28& \hfill & \hfill 54& \\ \hfill 8& \hfill 4& \hfill 36& & \hfill 70& \hfill \\ \hfill & \hfill & \hfill & & & \end{array}$$
(5.16)
Checks indicate that such a packing exists. It is known that the orthoplex bound cannot be achieved by 10 planes in $`G(3,2)`$, while the other cases are undecided. Our computer experiments strongly suggest that no set of 15 planes meets the simplex bound in $`G(5,2)`$.

On the other hand it is possible to find packings of 70 planes in $`G(8,4)`$ meeting the bound (5.14). With a considerable amount of effort (a long story, not told here), we determined several examples, of which the following is the most symmetrical. Let the coordinates be labeled $`\mathrm{},0,1,\mathrm{},6`$, and take two 4-planes generated by the vectors

$$\begin{array}{cccc}\{10000000,\hfill & 01000000,\hfill & 00100000,\hfill & 00001000\}\hfill \\ \{11000000,\hfill & 00101000,\hfill & 00010001,\hfill & 00000110\}\hfill \end{array},$$
(5.17)
respectively. We obtain 70 planes from these by negating any even number of coordinates, and/or applying the permutations $`(0123456)`$, $`(\mathrm{}0)(16)(23)(45)`$ and $`(124)(365)`$. The principal angles are $`0,0,\frac{\pi }{2},\frac{\pi }{2}`$; $`\frac{\pi }{4},\frac{\pi }{4},\frac{\pi }{4},\frac{\pi }{4}`$; or $`\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2}`$, so $`d_c^2=2`$, $`d_g^2=\pi ^2/4`$ (this is not even a local optimum with respect to geodesic distance). The full group, which is transitive on the planes, is generated by the above operations and by the $`8\times 8`$ Hadamard matrix

$$\frac{1}{\sqrt{8}}\left[\begin{array}{c}\\ ++++\\ ++++\\ ++++\\ ++++\\ ++++\\ ++++\\ ++++\end{array}\right],\begin{array}{c}\text{(}\text{ means }1\\ +\text{ means }+1\text{)}\end{array}.$$

It has shape $`2^8𝒜_8`$ and order 5160960 ($`𝒜_n`$ denotes an alternating group of degree $`n`$).

A set of 28 planes in $`G(7,3)`$ meeting the simplex bound can be obtained as follows. Label the coordinates $`0,\mathrm{},6`$, and let $`v_r^\pm `$ have components 1 at $`2^r(mod7)`$, $`\pm \sqrt{2}`$ at $`3.2^r(mod7)`$, and 0’s elsewhere, for $`r=0,1,2`$. Then the vectors $`v_0^\pm `$, $`v_1^\pm `$, $`v_2^\pm `$ (product of signs is even) span four planes, and the full set of 28 is found by cycling the seven coordinates.444Note added April 1996. P. W. Shor and the third author have recently discovered that the packing of 70 planes in $`G(8,4)`$ can be generalized to a packing of $`m^2+m2`$ planes in $`G(m,m/2)`$ meeting the orthoplex bound, whenever $`m2`$ is a power of 2; and that the packing of 28 planes in $`G(7,3)`$ can be generalized to a packing of $`p(p+1)/2`$ planes in $`G(p,(p1)/2)`$ meeting the simplex bound, whenever $`p`$ is a prime which is either 3 or congruent to $`1`$ modulo 8. These constructions will be described elsewhere.

On the other hand, in spite of much effort, we have not been able to find a set of 40 planes in $`G(6,3)`$ that meet the bound. If this were possible, we would obtain $`d_c^2=1.5`$, whereas the best we have been able to achieve is 1.49977. We can get 1.5 with $`N=34`$ planes, and 1.49977 with $`35N40`$, but only 1.4297 with 41 planes, which suggests that 1.5 might be possible with 40 planes. The discussion given earlier shows that this question could in principle be settled by seeing if a certain set of quadratic equations with integer coefficients has a real solution. One formulation leads to 500 quartic equations in 360 unknowns, so this approach seems hopeless at the present time.

It is not difficult to invent ways to construct sets of subspaces. Several examples are given below. However, many promising ideas have proved useless when confronted with the results found by our program. For this reason we have included two extensive tables of the best chordal distances that we have found. Table 3 gives $`d_c^2`$ for packings in $`G(m,2)`$, with $`m10`$. This will serve as a standard against which readers can test their own constructions.

$$\begin{array}{c}\\ \text{Table 3 about here}\end{array}$$

#### General constructions.

The following are some promising general constructions.

(i) A skew-symmetric conference matrix555For a list of the known orders of such matrices, see Table 7.1 of . of order $`4a`$ yields a set of $`4a`$ unit vectors in $`C^{2a}`$ with Hermitean inner products $`\pm i/\sqrt{4a1}`$ (, Example 5.8), hence a set of $`4a`$ planes in $`G(4a,2)`$ with $`d_c^2=4(2a1)/(4a1)`$ that meet the simplex bound.

(ii) Use the planes defined by the $`n`$-faces of a regular $`m`$-dimensional polytope, or the Voronoi or Delaunay cells of a lattice, etc. The initial results from this idea have been disappointing. For example, the 96 two-dimensional faces of the 24-cell in $`R^4`$ define 16 different planes, forming a packing with $`d_c^2=8/9`$, inferior to the best such packing (cf. Table 3).

(iii) Use the minimal vectors in a complex or quaternionic lattice to obtain packings in $`G(2a,2)`$ or $`G(4a,4)`$. For example, the 54 minimal vectors of the lattice $`E_6^{}`$, regarded as a three-dimensional lattice over the Eisenstein integers (, p. 127), produce a packing of nine planes in $`G(6,2)`$ that meets the simplex bound.

(iv) Restrict the search to generator matrices of 0’s and 1’s (as in (5.17)) or $`+1`$’ and $`1`$’s, or even to rows that are blocks in some combinatorial design, or vectors in some error-correcting code (see the example in Section 6).

(v) Let $`C`$ be an error-correcting code of length $`m`$ over $`GF(2^n)`$, for example a Reed-Solomon code (cf. ). If we regard $`GF(2^n)`$ as a vector space of dimension $`n`$ over $`GF(2)`$, each codeword yields an $`n\times m`$ matrix whose elements we may take to be $`+1`$’s and $`1`$’s. After discarding those of rank less than $`n`$, and weeding out duplicates, we obtain a packing in $`G(m,n)`$. The hexacode (, p. 82), for example, a code of length 6 over $`GF(4)`$ containing 64 codewords, produces 28 distinct planes in $`G(6,2)`$ with $`d_c^2=3/4`$. Unfortunately Table 4 shows that the record is 1.2973.

(vi) Choose a group $`G`$ with an $`m`$-dimensional representation, and a subgroup $`H`$ of index $`N`$ with an $`n`$-dimensional irreducible representation. Find an $`n`$-dimensional subspace $`VR^m`$ invariant under $`H`$, and take its orbit under $`G`$. Many of the conference matrix constructions of line-packings described in Section 6 are of this type (using $`G=L_2(q)`$).

In some cases we have also used the optimizer to search for packings with a specified group. The following is a packing of 28 planes in $`G(8,2)`$ meeting the simplex bound that is an abstraction of a configuration found by the computer when searching for packings in $`R^8`$ invariant under the permutation $`(0)(1,2,3,4,5,6,7)`$. Let $`R:C^4R^8`$ map $`(v_1,v_2,v_3,v_4)`$ to $`(\mathrm{Re}v_1,\mathrm{Im}v_1,\mathrm{Re}v_2,\mathrm{},\mathrm{Im}v_4)`$, and let $`\alpha =e^{2\pi i/7}`$. The 28 planes are spanned by the following pairs of vectors:

$$\begin{array}{ccc}R(0,\alpha ^k,\alpha ^{2k},\alpha ^{4k}),\hfill & R(0,i\alpha ^k,i\alpha ^{2k},i\alpha ^{4k}),\hfill & 0k6,\hfill \\ R(1,0,\alpha ^{2k},\alpha ^{4k}),\hfill & R(i,0,i\alpha ^{2k},i\alpha ^{4k}),\hfill & 0k6,\hfill \\ R(1,\alpha ^k,0,\alpha ^{4k}),\hfill & R(i,i\alpha ^k,0,i\alpha ^{4k}),\hfill & 0k6,\hfill \\ R(1,\alpha ^k,\alpha ^{2k},0),\hfill & R(i,i\alpha ^k,i\alpha ^{2k},0),\hfill & 0k6.\hfill \end{array}$$

(The pattern of zeros and signs here suggests the Tetracode , p. 81.)

#### An application: apportioning randomness.

The following is a potential application of packings in $`G(m,n)`$ for use in apportioning randomness, in the sense of producing large numbers of approximately random points from a few genuinely random numbers. We illustrate using the example of 70 planes in $`G(8,4)`$ described above. Let $`A_1,\mathrm{},A_{70}`$ be generator matrices for them. Suppose an exclusive resort wishes to distribute random points in $`R^4`$ to its 70 guests, perhaps for use as garage door keys. Let $`x=(x_1,\mathrm{},x_8)`$ be a vector of eight independent Gaussian random variates. Then the hotel would assign $`y_i=A_ix^{tr}`$ to its $`i`$-th guest. By maximizing the chordal distance between the planes we have minimized the correlation between the $`y_i`$.

## 6 Packing lines in higher dimensions

Table 4 shows the maximal angular separation found for packings of $`N50`$ lines in $`G(m,1)`$ for $`m9`$.

The following packings of $`N`$ lines in $`G(m,1)`$ for $`m10`$ achieve either the simplex or orthoplex bounds: $`m=3:N=3,4,6,7`$; $`m=4:N=4,5,11,12`$; $`m=5:N=5,6,10,16`$; $`m=6:N=6,7,16,22`$; $`m=7:N=7,8,14,28`$; $`m=8:N=8,9`$; $`m=9:N=9,10,18,4648`$, $`m=10:N=10,11,16`$. Most of these are described below.

As mentioned in Section 3, the best packing of 5 lines in $`G(3,1)`$ is a subset of the best packing of 6 lines. Table 4 shows similar phenomena in higher dimensions. For example, the putatively best packing of 48 lines in $`G(9,1)`$ is so good that we cannot do better even when up through 8 lines are omitted from it.

In the rest of this section we discuss some of the entries in Table 4 that have the largest symmetry groups.

$$\begin{array}{c}\\ \text{Table 4 about here}\end{array}$$

#### Constructions from lattices.

Let $`P_1,\mathrm{},P_a`$ be a set of mutually touching spheres in a $`d`$-dimensional lattice packings. If there are $`2N`$ further spheres in the packing, each of which touches all of $`P_1,\mathrm{},P_a`$, their centers form a $`(da+1)`$-dimensional antipodal spherical code with angular separation $`\mathrm{arcsec}(a+1)`$ (, p. 340, Theorem 1).

In particular, the entries labeled “1” in Table 4 are obtained from the centers of spheres that touch one sphere in the lattices $`D_4`$, $`D_5`$, $`E_6`$ and $`E_7`$. The entries labeled “2” are obtained from spheres that touch two spheres in the lattices $`D_4`$, $`D_5`$, $`E_6`$, $`E_7`$, $`E_8`$, $`\mathrm{\Lambda }_9`$ and the nonlattice packing $`P_{10b}`$ (cf. , Chap. 1, Table 1.2). We remark the set of 28 lines in $`G(7,1)`$ with angle $`\mathrm{arccos}(1/3)`$ obtained in this way from $`E_8`$ is known to be unique: see Chapter 14, Theorem 12 of . This configuration of lines is derived in a different way in .

It is interesting to compare Table 4 with the table of maximal sets of equiangular lines given in Lemmens and Seidel  and Seidel . Some arrangements appear in both tables, for example the set of 28 lines in $`G(7,1)`$ just mentioned. On the other hand the difference between the tables can be seen in dimension 8. The maximal set of equiangular lines that exists in $`G(8,1)`$ has size 28, with angle $`\mathrm{arcsec}3`$ (, Theorem 4.6). However, Table 4 shows that there is a set of 28 lines (not equiangular) in $`G(8,1)`$ with minimal angle arcsec 3.000511… . In view of the Lemmons-Seidel result (which uses the fact that in an $`m`$-dimensional equiangular arrangement of $`N`$ lines with $`N2m`$ the secant must be an odd integer) our set of 28 lines cannot be perturbed to give an equiangular set without decreasing the minimal angle.

#### From conference matrices.

It follows from Theorem 6.3 of  that if a symmetric conference matrix of order $`q+12`$ $`(mod4)`$ exists then there is an arrangement of $`q+1`$ equiangular lines in $`R^m`$, $`m=(q+1)/2`$, with $`d_c^2=(q1)/q`$, meeting the simplex bound. The corresponding entries are labeled “3” in the table. Our program was able to find these packings for every prime power $`q`$ of this form below 100 except for 49 and 81.

#### From diplo-simplices.

The entries labeled “4” are obtained by using all the vectors of shape $`\pm c(n^1,(1)^n)`$, where $`c=1/\sqrt{n(n+1)}`$. In the notation of , Chap. 4, these are the minimal vectors in the translates of the root lattice $`A_n`$ by the glue vectors $`[1]`$ and $`[1]`$. They form the vertices of a diplo-simplex , and have automorphism group $`2\times 𝒜_{n+1}`$. These packings also meet the simplex bound.

#### From codes.

Let $`C`$ be a binary code of length $`m`$, size $`M`$ and minimal distance $`d`$ which is closed under complementation. Writing the codewords as vectors of $`\pm 1`$’s, we obtain a packing of $`M/2`$ lines in $`G(m,1)`$ with $`d_c^2=4d(md)/m^2`$. For example, a shortened Hamming code of length 10 with $`M=32`$, $`d=4`$ gives a packing of 16 lines in $`G(10,1)`$ meeting the simplex bound. This construction provides a rich supply of good packings. The Nordstrom-Robinson code , , for example, yields a packing of 128 lines in $`G(16,1)`$ with $`d_c^2=15/16`$.

#### 40 lines in $`G(4,1)`$.

The 80 antipodal points are $`(\xi ^{\mu +1/2},0)`$, $`(0,\xi ^{\nu +1/2})`$, $`(a\xi ^{2\mu },b\xi ^{2\nu +1})`$, $`(a\xi ^{2\mu +1},b\xi ^{2\nu })`$, $`(b\xi ^{2\mu },a\xi ^{2\nu })`$, $`(b\xi ^{2\mu +1},a\xi ^{2\nu +1})`$, where $`\xi =e^{2\pi i/8}`$, $`\mu ,\nu =0,\mathrm{},3`$, $`a=2^{1/4}`$, $`b=\sqrt{1a^2}`$. The group of these 80 points is the group $`\frac{1}{2}(𝒟_{16}+𝒟_{16})2`$ of order 256 (see ). Here $`𝒟_m`$ denotes a dihedral group of order $`m`$.

#### 22 lines in $`G(6,1)`$.

We use the vertices of a hemi-cube (an alternative way to describe the 16-line packing) together with the six coordinate axes.

#### 63 and 64 lines in $`G(7,1)`$.

The putatively best packing of 63 lines in $`G(7,1)`$ (just beyond the range of Table 4) has angular separation $`60^{}`$, and is formed from the 126 minimal vectors of the lattice $`E_7`$. The automorphism group of this set of 126 points is the Weyl group $`𝒲(E_7)`$, of order $`2^{10}.3^4\mathrm{.5.7}`$ (cf. , Chap. 4, §8.3).

The best packing found of 64 lines in $`G(7,1)`$ also has an unusually large group, of order $`2|𝒲(E_6)|=2^8.3^4.5`$. This packing can be obtained as follows. The largest gaps between the packing of 63 lines occur in the directions of the minimal vectors of the dual lattice $`E_7^{}`$. Let $`u`$ be such a vector, with $`u=\text{“sevenrm1}\text{/}\text{“sevenrm2}v_6`$, $`v_6E_7`$, $`v_6v_6=6`$. Taking a coordinate frame in which the last coordinate is in the $`v_6`$ direction, we obtain the minimal vectors of $`E_7`$ in the form $`(w,0)`$, $`wE_6`$, $`ww=2`$ (72 vectors), and $`(x,\pm \sqrt{2/3})`$, $`xE_6^{}`$, $`xx=4/3`$ (54 vectors). We now adjoin $`\pm \sqrt{6}`$ and rescale by multiplying the last coordinate by $`\sqrt{3/4}`$, obtaining our final configuration of vectors $`(0,\pm 3/\sqrt{2})`$, $`(w,0)`$, $`(x,\pm 1/\sqrt{2})`$, which is isomorphic to what the computer found. The rescaling has compressed the $`E_7`$ lattice in the $`v_6`$ direction until the angle between $`v_6`$ and the $`(x,1/\sqrt{2})`$ vectors is the same as the minimal angle between the $`(x,1/\sqrt{2})`$ and $`(w,0)`$ vectors. This angle, the minimal angle in the packing, is $`\mathrm{arccos}(\sqrt{3/11})=58.5178^{}`$.

#### 36 lines in $`G(8,1)`$.

The 72 antipodal points consist of all permutations of the two vectors $`c(7^2,2^7)`$ and $`c(7^2,2^6)`$, where $`c=1/\sqrt{126}`$. These are the minimal vectors in the translates of the $`A_8`$ lattice by the glue vectors $`[2]`$ and $`[2]`$.

The configuration of 16 lines in $`G(5,1)`$ is similarly obtained from the translates of $`A_5`$ by $`[1]`$, $`[3]`$ and $`[5]`$.

#### 39 lines in $`G(12,1)`$.

Take the 13 lines of a projective plane of order 3, described by vectors of four 1’s and nine 0’s, and replace two of the 1’s by $`1`$’s in all possible ways, obtaining 78 vectors in $`R^{12}`$ with angle arcsec 4. The group is $`2\times L_3(3)`$, of order $`2^53^313`$.

If instead we replace any odd number of the 1’s by $`1`$’s, we obtain a putatively optimal packing of 52 lines in $`G(13,1)`$, with angle arcsec 4 and group $`2\times L_3(3).2`$ of order $`2^63^313`$.

As in all these examples, we were amazed that the program was able to discover such beautiful configurations.

### Acknowledgements

We are grateful to Daniel Asimov, Andreas Buja, Dianne Cook, Simon Kochen, Warren Smith, and especially Colin Mallows for helpful discussions, and to David Applegate for help in solving the matching problem described in Section 4.


format: text

text:
# Reviewing Pulsar Spectra

## 1 Problems solved?

Pulsar research must be considered - 35 years after the detection of pulsars - a mature science, where the basic questions have been raised and discussed. One would hope that many if not all generic and important problems have found some kind of answer and that scientific work can concentrate now on specific details requiring more in depth investigation. We know, however, that this picture is not true. Even well studied areas did in the past not always lead to a general accepted model and some were investigated at the beginning with enthusiasm but left behind.

This paper will concentrate on one narrow topic, pulsar radio spectra. It is the attempt to work out, what features are now generally accepted, but also what features are still in discussion after so many years of pulsar research. In this sense the paper will enlighten in a personal view those aspects which are still under discussion.

## 2 The spectrum, indicator of the emission mechanism

The spectrum, i.e. the total emitted energy in specific radio bands, was considered from the beginning as an easy to measure generic characteristic of the emission mechanism, especially since one could hope to learn something about an exciting new type of coherent amplification. Comella (1971) collected and published the first data points from measurements made at Arecibo, a work which was carried on by Sieber (1973) covering a wider compilation of observations from different observatories. Later on comprehensive data collections were published by Izvekova et al. (1981) and Lorimer et al. (1995) who extended the compilation to several hundreds of pulsars, albeit concentrated on few frequency bands. More recently wide collections of data were published by Malofeev (2000) and Maron et al. (2000).

From the beginning it was clear that pulsar spectra are characterizes by a ”steep” spectral index, steeper than e.g. measured for supernova remnants, which explains why they were detected at rather low - below 1 GHz - frequencies. There were in addition indications that the spectra even steepen in some cases at higher frequencies beyond 1 GHz, although the measurement uncertainty was high. Undoubtedly, some pulsars showed a ”low-frequency turnover”, a decrease of the emitted energy at low frequencies, typically around 100 MHz, whereas other pulsars had a straight spectrum down to the lowest observable frequencies. The cause of this decrease - some kind of absorption in the magnetosphere or a loss of efficiency of the emission mechanism was discussed - remained unclear.

The problem with all the measurements was a high variability of the pulsar emission on many time scales, which resulted in a high uncertainty of the data points and which made it questionable, if a constant, time invariable spectrum exists at all. Investigations of pulsar spectra have, therefore, to face this problem and to look into the various possible contributions to the variability.

## 3 Variability

### 3.1 Intrinsic variability

Emission within a pulse during the pulse window is - depending on the pulsar - highly structured showing sometimes features as short as microseconds (microstructure) or even nanoseconds. Most pulsars have clearly visible sub-pulses extending over milliseconds, which in some rare cases are well organized with a shift from pulse to pulse, the so called drifting sub-pulses, Fig. 1.

This effect can now be quantitatively explained by the occurrence of an ExB force in the magnetosphere working on the out streaming particles. The drift bands cause a well defined periodic modulation (variability) of the received emission from pulse to pulse called P3 with characteristic time scale of several periods. Interestingly, some pulsars exhibit periodic modulations without recognizable drifting sub-pulses, where the duration extends over many periods, a kind of memory of the radiation mechanism which seems to be poorly understood so far.

Also poorly understood seems to be the producing mechanism for a sudden drop of the emission to zero (at least below the detection limit) called nulling, an effect which can be found for many pulsars. These nulls can stretch over many periods, even dominating the emission pattern, so that the emission occurs only in well defined bursts; an effect, which obviously results in a strong modulation of the radiation which may extend over many minutes of time.

Measurements show that the mentioned intensity variations are essentially broadband ($`>100\mathrm{M}\mathrm{H}\mathrm{z}`$) and with high probability related to the emission mechanism or the structure of the magnetosphere or the surface of the neutron star. The time scales extend from nanoseconds to minutes. It is the hope (we will come to this point later) that averages of the received signal over many minutes will ”smooth out” all these variations.

There is less confidence with respect to another effect, which was called mode switching: In contrast to the well known behavior that the average of many pulses results in a stable average pulse profile, some pulsars show sudden variations of the average profile normally correlated with a change in intensity (Fig. 2). The times of variation are unpredictable and may occur suddenly even after e.g. several hours (Fig. 3).

For mode switching pulsars, time stable intensity values are obtainable only after some hours of integration, depending on the degree of accuracy which one would like to reach (see e.g. the investigation of Helfand et al. (1975) on the stability of integrated pulse profiles). Nevertheless, the hope is, that the intrinsic intensity variations can be averaged out by measurements lasting long enough - at least many minutes, for some pulsars even several hours. Least stable seem to be the mode switching pulsars whose long-term behavior has not yet been investigated very well.

### 3.2  External influence

#### 3.2.1  Diffractive interstellar scintillation DISS

Scheuer (1968) showed as early as 1968 that interstellar scintillation may be the cause of the modulation of the received pulsar signal on time scales around several minutes, mainly because the small extent of the emission region makes pulsars ideal ”point sources”. Rickett (1969) could later proof that indeed features of strong and weak scintillation can be found in the received emission of pulsars, a decorrelation of the received signal in frequency and time being the clearest signature of this scintillation process, now being called diffractive interstellar scintillation DISS. The decorrelation is best represented by way of ”dynamic spectra” (Fig. 4), where the average emission of one or several pulses is analyzed in a spectrometer and displayed as intensity contours. Successive spectra then depict the development of the structure in time.

Typical time scales for diffractive scintillation range from seconds to hours and frequency scales from about 100 kHz to a few 100 MHz. Since scattering parameters are now known for many pulsars, at least in the strong scintillation regime - albeit by far not for all and not for each frequency range - it is in principle possible to remove the strong intensity variations caused by DISS. Carefully planned observations, extending over much more than the decorrelation bandwidth, and/or with durations much longer than the decorrelation time should ”average out” the variations so that the intrinsic emitted signal becomes measurable.

Alternatively one could take many observations of the emitted energy at a certain frequency - even if bandwidth and duration are smaller than the decorrelation scales - and average these scintillation modulated signals afterwards. The averaging process should smooth out the variations if enough samples are taken. Interestingly enough, this was more or less the procedure which was taken intuitively by the first observers, who just collected the measurements of many authors or observing periods.

#### 3.2.2  Refractive interstellar scintillation RISS

For years it was known that pulsars show intensity variations even if diffractive scintillation effects are removed. The cause of the remaining long-term variations was unclear until Sieber (1982) could show, that they are correlated with the dispersion measure and due to the interstellar medium. The theoretical foundation was laid by Rickett, Coles and Bourgois in 1984, who showed that focussing and defocussing by features large compared to the Fresnel scale are the origin, an effect which is now called refractive interstellar scintillation RISS. The effect is rather broadband and the times scales rather large, extending perhaps even up to months or years as demonstrated e.g. by long-term observations of Stinebring et al. (2000) and Ramesh Bhat et al. (1999).

The results of Ramesh Bhat et al. (1999) shook, unfortunately, the confidence that diffractive scintillation effects may be easily removable by averaging over the decorrelation time or frequency scale, because they showed that the time and frequency scale is also variable, being influenced by the longer lasting RISS. Only really long-term observations are apparently able to tell the truth about the emitted signal strength.

The measurements of Stinebring et al. (2000) contained also good news: Their observations demonstrate quite clear that some pulsars, like PSR B0736-40, B1641-45 and B1859+03 exhibit very stable intensity over many years of observations. This means that the intrinsic emission at least for these pulsars must be stable. The observations give a clear hint, that the intrinsic emission of perhaps all pulsars is stable on long terms (intrinsic effects being excepted as demonstrated in Sec. 3.1). So that we have good hope that we can talk about stable intrinsic (emitted) pulsar radio spectra (on time scales as discussed above); which means that the original, naive, approach to the problem of pulsar spectra - being based on the assumption of stable intrinsic spectral behavior - has found a justification. There remains the problem, that very few pulsars are known so far whose emission has been measured to be stable on long terms. Far more long-term observations are needed to create a sound basis for accurate pulsar spectra determinations.

## 4  Spectra

### 4.1  Spectral index

0n the basis of the discussion in the preceeding sections, it seems appropriate to talk about pulsar spectra as a well-defined intrinsic characteristic of pulsars. Observations of the first years showed already the steep character of the spectra (Sieber 1973), with an exponent of -1.62 for the high-frequency, power law part of the spectrum (not taking into account the low-frequency turnover found for some pulsars and without considering the high-frequency steepening). These first investigations were based on few sources (27), a number, which was increased later on by Izvekova et al. (1981) (52 sources) and then in a big step to 280 by Lorimer et al. (1995), Malofeev (1996, 2000) and Maron et el. (2000).

The computed mean spectral indices vary for these investigations between -1.6 and -2.0, underlining the steep character of the emitted spectra, but showing also a remaining uncertainty in the determination of the average value (which makes it not very sensible to give error bars). There is a wide scatter of measured spectral indices apparent for individual pulsars, ranging from 0.0 (or even slightly positive) to -3.1, which may explain also the scatter in the compiled average spectral index. There is also some uncertainty due to the fact that the frequency range used for the determination of the exponent is often limited either by lack of data or because one has to exclude either the low-frequency or the high–frequency turnover.

The spectral index is only weakly correlated with any other measurable pulsar parameter, as e.g. P, dP/dt etc., the clearest correlation being given to the characteristic age (Lorimer et al. 1995). It remains to theories of the emission mechanism to explain the physical reasons for the differences in spectral index for different pulsars (See however section 6).

### 4.2  Low-frequency turnover

A high percentage of those (normal) pulsars, where flux values below 100 MHz have been measured, show a low-frequency drop (turnover) in intensity (see Fig. 5); and it is known that correlations between for example the maximum frequency and the pulsar period exist.

Several mechanism have been proposed to explain this drop in intensity, often related to the radiation mechanism. Since the debate on the most probable emission mechanism is still ongoing, there is also still no definite answer to the cause of this low-frequency turnover.

It is interesting to note, that some pulsars on the contrary show definitely no low-frequency turnover, although accurate measurements are available for them down to very low frequencies (10 MHz). Famous examples are the Crab pulsar (Fig. 6) and the millisecond pulsar PSR B1937+214 (Erickson & Mohany 1985).

Explanations for the low-frequency turnover should also clearly state, why some pulsars obviously show increasing intensity down to very low frequencies.

### 4.3  Break in the power law spectrum

There were already at the beginning of the investigation of pulsar spectra indications, that also the high-frequency part of the spectrum may not completely be describable by just one power law. Very few pulsars showed indications of a steepening of the spectrum (see e.g. Fig. 7), which could be fitted by two power laws. The fact that pulsar emission is highly variable on many time scales (Sec. 3), makes it difficult to firmly proof such a break in the spectrum and to separate clearly intrinsic effects from e.g. long-time refraction effects. It should, however, be mentioned, that a break in the spectrum is a natural consequence of some theoretical radiation emission mechanisms; see e.g. Ochelkov & Usov (1984).

Most interesting appears also the fact that some pulsars show obviously an increase of their emission intensity at very high frequencies; this effect will be discussed in detail in the contribution of Wielebinski (2002) in these proceedings.

## 5  Millisecond pulsars

The detection of millisecond pulsars opened up a new era in the investigation of pulsars. Millisecond pulsars occupy a different region in the parameter space and build therefore a severe test for every theory of radio emission. Comparisons between the characteristics of normal pulsars and millisecond pulsars were therefore eagerly expected. Kramer et al. (1998) published the spectra of 32 millisecond pulsars and found an exponent of = -1.8 if no distance limit was applied. Nearer pulsars (distance less then 1.5 kpc) showed a bit flatter spectrum with = - 1.6. They compared this result with 346 normal pulsars, which showed an exponent of = -1.6 for all distances and of = - 1.7 for distances $`<1.5\mathrm{kpc}`$. In correspondence to the result of Toscano et al. (1998), who measured an index of = -1.9 for southern millisecond pulsars, one comes to the conclusion that there exists no real difference in the spectral index between normal pulsars and millisecond pulsars. It seems that the same radiation mechanism works in both environments.

### 5.1  Low-frequency turnover

Kuzmin & Losovsky (2001) were able to extend the spectral frequency range of quite a collection of millisecond pulsars down to 100 MHz. Interesting enough no low-frequency turnover (drop at low frequencies) was detected, with one exception, PSR J1012+5307. Together with the fact that the millisecond pulsar PSR B1937+214 shows a straight spectrum down to 10 MHz, there is the indication, that low-frequency turnovers might be rarer in millisecond pulsars than in normal pulsars.

## 6  Geometry

The geometry of the hollow-cone must have a direct impact on the spectrum, if the radio emission is - as we believe - concentrated in a hollow cone with a frequency dependent width. Kramer et al. (1994) and Sieber (1997) could demonstrate that the observed spectrum may be indeed steeper than the emitted one, the steepening being caused by the shrinking of the hollow-cone width at high frequencies.

An emitted spectrum with exponent = -1.5 could - depending on the specific geometry - be observed as a spectrum with = -2. This would mean, that the radiation mechanism itself must not produce such steep spectra as we observe them.

In addition, the received intensity may depend on several more geometrical factors, e.g. the curvature radius of the field lines, the emission hight over the neutron star surface and the orientation of magnetic field axis and rotation axis. Kuzmin & Solov’ev (1986) worked out an curvature radiation emission model which reflects quite well as an example the spectral characteristics of PSR B1133+16, including the low-frequency turnover and taking into account all these geometrical effects.

It appears obvious that geometrical effects influence the shape of the received spectrum; but it is much harder to determine the exact impact in individual cases, since then the precise geometry has to be known, e.g. the angle between rotation axis and magnetic field , as well as the impact angle . Normally the geometry is derived from measurements of the polarization angle sweep across the profile; these measurements are, however, often spoiled by orthogonal polarization angle jumps and may give unreliable results. Most promising are polarization angle histograms which often allow to determine at least the angle sweep with good precision.

## 7  Conclusions

From the arguments presented one is led to the conclusion that astonishingly few problems found an answer so far. One of these is obviously the question of the stability of the intrinsic spectrum which seems to be settled in the sense, that we can now indeed assume a stable emitted spectrum being heavily modulated by transmission effects on all time scales extending up to months. Also the intrinsic variations manifested as drifting sub-pulses and P3 modulation can now be explained more and more in the framework of geometrical mapping to the polar cap of the neutron star.

To the still debated problems belongs e.g. the question of the cause of the low-frequency turnover in the spectra of pulsars - normal and millisecond - the cause of the high-frequency break in some spectra and the turn-up at extreme high frequencies. And finally it would be of great help if the geometry of the emission geometry of more and more pulsars could be determined with higher confidence.

Acknowledgments: The support of the Wilhelm und Else Heraeus-Stiftung is very much acknowledged. This work was also supported by INTAS grant no. 2000-849.


text:
# Coreless vortex ground state of the rotating spinor condensate

## I Introduction

The quantization of circulation for a Bose-Einstein condensate (BEC) is well known . This quantization implies a crucial role for quantized vortices as carriers of angular momenta. If the condensate rotates with angular frequency $`\mathrm{\Omega }`$, the ground state of the system is a configuration of quantized vortices and the number of vortices increases with the rotation frequency . Energetically a state with a single vortex becomes favored only after $`\mathrm{\Omega }`$ exceeds certain threshold $`\mathrm{\Omega }_c`$ that depends on the system parameters in a fairly complicated manner .

While these properties are well studied, the creation of spinor condensates  has added an important extra twist. In a spinor (or vectorial) condensate with spin $`S`$ the spin-degree of freedom is not frozen and the condensate order parameter is a vector $`\zeta `$ with $`2S+1`$ components. This freedom with the spin-degree has some intriguing consequences. For example, in inhomogeneous magnetic field this vectorial nature of the order parameter can even imply vortex ground states . Crucially for our purposes, such circulation of an order parameter does not have to be quantized . If circulation is not quantized angular momenta of the system does not have to be carried by quantized vortices.

In this paper we address the question of the ground state of the rotating spinor condensate. We show that the ground state of the ferromagnetic spinor condensate at small rotation frequencies is a coreless (non-singular) vortex and we analyze its properties. Coreless vortices have also been studied by Mizushima et al.  who investigated rotating spinor condensates with a fixed magnetization.We mainly concentrate on the pure rotational ground state with no external magnetic field. Furthermore, we study the rotating spinor condensate with polar (anti-ferromagnetic) interactions and show how significantly the behavior of the system changes when only the sign of the antisymmetric interaction parameter is altered. In both cases there is a critical rotation frequency above which the ground state carries non-zero angular momentum. However, in the ferromagnetic case the amount of angular momentum can be infinitesimal whereas the corresponding polar condensate has a quantized angular momentum (in trap units) $`1/2`$ per particle. The vortex structure of a polar spinor BEC has been previously studied by Yip  focusing for a single constant rotation frequency.

## II Model

The Hamiltonian for the three-dimensional spin-1 condensate in a frame rotating with frequency $`\mathrm{\Omega }`$ is

$`H`$ $`={\displaystyle }d^3𝐫({\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_i+{\displaystyle \frac{m}{2}}(\omega _{tr}^2r^2+\omega _z^2z^2)\mathrm{\Psi }_i^{}\mathrm{\Psi }_i`$ (3)
$`\mathrm{\Omega }\mathrm{\Psi }_i^{}\widehat{L}\mathrm{\Psi }_i+{\displaystyle \frac{\lambda _s}{2}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}\mathrm{\Psi }_j\mathrm{\Psi }_i`$
$`+{\displaystyle \frac{\lambda _a}{2}}{\displaystyle \underset{\alpha }{}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}(F_\alpha )_{ik}(F_\alpha )_{jk}\mathrm{\Psi }_k\mathrm{\Psi }_l),`$
where $`m`$ is the atomic mass, $`\omega _{tr}`$ and $`\omega _z`$ are the trapping frequencies, $`\widehat{L}`$ is the angular momentum operator, and $`𝐅_\alpha `$ is the angular momentum matrix. The parameters $`\lambda _s`$ and $`\lambda _a`$ are system specific parameters that depend on the two different scattering lengths $`a_2`$ and $`a_0`$: $`\lambda _s=\frac{4\pi \mathrm{}^2(a_0+2a_2)}{3m}`$ and $`\lambda _a=\frac{4\pi \mathrm{}^2(a_2a_0)}{3m}`$. In the mean field limit our aim is to find the fields $`\mathrm{\Psi }_i(𝐫)=\psi _i(𝐫)`$ that minimize this Hamiltonian. Minimization leads to the three-component Gross-Pitaevskii equations . We choose $`\omega _z=(2\pi )1480`$ Hz and $`\omega _{tr}=(2\pi )62`$ Hz. When the trap is filled with $`1.7\times 10^4`$ Rb87 atoms the condition $`\mu <\mathrm{}\omega _z`$ is satisfied. Consequently, in $`z`$-direction we assume the system to be in the ground state of a parabolic trap with no dynamics and our numerical calculation becomes effectively two-dimensional. This assumption scales the scattering lengths by the constant $`\eta =\sqrt{\frac{\omega _z}{2\pi \omega _{tr}}}`$. For convenience we choose the unit of length as $`L_0=\sqrt{\frac{\mathrm{}}{m\omega _{tr}}}`$ and the unit of time as $`\tau _0=1/\omega _{tr}`$.

It is well known that the ground state structure $`\psi =(\psi _1,\psi _0,\psi _1)^T=\sqrt{n}\zeta `$ depends on the sign of the $`\lambda _a`$ . In the absence of a magnetic field, if $`\lambda _a>0`$ the energy is minimized (with a suitable choice of axis) with a non-magnetized spinor $`\zeta =(0,1,0)^T`$. This state is referred to as the polar state. If $`\lambda _a<0`$, a fully magnetized ferromagnetic state $`\zeta =(1,0,0)^T`$ is favored. In the absence of external magnetic field we can freely rotate the order parameter (globally) without any physical changes. In particular we can rotate the ferromagnetic ground state and get an identical order parameter expressed in terms of the Euler angles $`(\alpha ,\beta ,\tau )`$ of the rotation

$`\psi =\sqrt{n}e^{i\tau }\left(\begin{array}{c}e^{i\alpha }\mathrm{cos}^2\beta /2\\ \sqrt{2}\mathrm{cos}\beta /2\mathrm{sin}\beta /2\\ e^{i\alpha }\mathrm{sin}^2\beta /2\end{array}\right).`$ (7)
The superfluid velocity is defined by $`𝐯_𝐬=i\frac{\mathrm{}}{m}\zeta ^{}\zeta `$. If we make the rotation local and choose $`\tau =\alpha =\varphi `$, with $`\varphi =\mathrm{tan}^1(x/y)`$, the velocity field of the condensate is given by $`𝐯_𝐬=\frac{1}{r}(1\mathrm{cos}\beta )\widehat{\varphi }`$ . If $`\beta (r)=0`$ the velocity field vanishes, but as $`\beta (r)`$ increases the velocity field starts to resemble that of an ordinary singular vortex. Such a structure is called a coreless vortex which was first predicted by two studies for superfluid 3He. For the Mermin-Ho (MH)  vortex the bending angle $`\beta `$ must be $`\frac{\pi }{2}`$ at the boundary of the condensate and for the Anderson-Toulouse (AT)  vortex $`\beta `$ must be $`\pi `$. For the MH vortex in liquid 3He the boundary value is imposed by the walls of the container .

## III Ferromagnetic ground state

The angular momentum per particle in a condensate with a coreless vortex can be anything between $`0`$ and $`2`$. Also, a coreless vortex can be continuously transformed into a state with vanishing angular momentum. This happens by simply transforming the function $`\beta (r)`$ into zero. Due to these reasons it seems plausible that when $`\lambda _a<0`$, a coreless vortex might play a role as a carrier of angular momenta at small values of $`\mathrm{\Omega }`$. We have verified this reasonable conjecture by solving the relevant two-dimensional GP equations numerically for $`1.7\times 10^4`$ Rb87 atoms without any additional assumptions:

$`i\mathrm{}{\displaystyle \frac{\psi _1}{t}}`$ $`=`$ $`\psi _1+\lambda _a(\psi _0^2\psi _1^{}+|\psi _1|^2\psi _1+|\psi _0|^2\psi _1`$ (9)
$`|\psi _1|^2\psi _1)`$
$`i\mathrm{}{\displaystyle \frac{\psi _0}{t}}`$ $`=`$ $`\psi _0+\lambda _a\left(2\psi _1\psi _1\psi _0^{}+|\psi _1|^2\psi _0+|\psi _1|^2\psi _0\right)`$ (10)
$`i\mathrm{}{\displaystyle \frac{\psi _1}{t}}`$ $`=`$ $`\psi _1+\lambda _a(\psi _0^2\psi _1^{}+|\psi _1|^2\psi _1+|\psi _0|^2\psi _1`$ (12)
$`|\psi _1|^2\psi _1),`$
where $`=\frac{\mathrm{}^2}{2m}^2+V_{trap}\mathrm{\Omega }\widehat{L}+\lambda _s\left(|\psi _1|^2+|\psi _0|^2+|\psi _1|^2\right)`$. The lowest energy configuration at small values of $`\mathrm{\Omega }`$ is a coreless vortex (see an example of the superfluid velocity field in Fig. 1). A priori, there is no reason to expect a critical rotation frequency for the existence of a coreless vortex. Nevertheless, by solving the GP equations we saw that a critical frequency $`\stackrel{~}{\mathrm{\Omega }}_c0.06`$ had to be exceeded before a finite sized coreless vortex became energetically favored. Below this critical frequency energy was minimized with the function $`\beta (r)`$ being zero everywhere.

Assuming the existence of a coreless vortex the Hamiltonian in trap units is

$`H`$ $`{}_{CV}{}^{}[n,\beta ]={\displaystyle }d^2rn[{\displaystyle \frac{1}{2}}r^2+{\displaystyle \frac{\lambda _s+\lambda _a}{2}}n+{\displaystyle \frac{1}{2n}}|\sqrt{n}|^2`$ (14)
$`+{\displaystyle \frac{1}{4}}|\beta |^2+{\displaystyle \frac{3+\mathrm{cos}^2\beta 4\mathrm{cos}\beta }{4r^2}}2\stackrel{~}{\mathrm{\Omega }}\mathrm{sin}^2\beta /2],`$
where $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }/\omega _{tr}`$. When the rotation frequency is low and the density is large, it is reasonable to assume a density distribution almost independent of the function $`\beta `$ and that this density distribution has the familiar inverted parabola shape of the Thomas-Fermi approximation . The $`\beta (r)`$ minimizing the energy functional (14) must satisfy the following differential equation

$`{\displaystyle \frac{^2\beta }{r^2}}`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{\beta }{r}}{\displaystyle \frac{\frac{\beta }{r}\frac{n}{r}}{n}}+{\displaystyle \frac{4\mathrm{sin}\beta \mathrm{sin}2\beta }{2r^2}}`$ (17)
$`2\stackrel{~}{\mathrm{\Omega }}\mathrm{sin}\beta .`$
This equation is solved numerically to get the lowest energy configuration. Some results are shown in Figure 2. By using the same parameters as in the numerical simulations the energetical studies give $`\stackrel{~}{\mathrm{\Omega }}_c=0.058`$ for the critical rotating frequency. This result is extremely close to the corresponding value from the GP equations. Also, the values of angular momenta as a function of $`\stackrel{~}{\mathrm{\Omega }}`$ are nearly identical as can be seen in Figure 3.

To gain more insight, Eq. (14) can also be studied variationally. With a trial $`\beta (r)=r/L`$ the resulting integrals can be solved. While the resulting expression is quite complicated, it is remarkable that such a simple variational model gives results that are qualitatively the same as the results based on the Eq. (17). It turns out that the the existence of an energy minimum depends (in trap units) on the magnitude of the number $`A=R_{TF}^2\stackrel{~}{\mathrm{\Omega }}`$. Energy mimimum only exists if $`A>3.05\pm 0.05`$. Thus the critical frequency predicted by the variational model with the parameters used before is $`\mathrm{\Omega }_c=0.064\pm 0.001`$, i.e. very close to the value calculated with more accurate methods. As the critical frequency is inversely proportional to the square of the Thomas-Fermi radius we see that the non-zero critical rotation frequency is a finite size effect. This property is also confirmed by the numerical solutions of the GP equations.

### A Ground state as a function of $`\stackrel{~}{\mathrm{\Omega }}`$

As the rotation frequency is increased, the slope of $`\beta `$ increases initially (to a good accuracy) linearly with the rotation frequency. This can be easily understood since initially a larger slope implies larger populations for the state with the highest winding number, namely $`m=1`$ state. But if the slope becomes too large $`\beta (r)`$ will be larger than $`\pi `$ inside the condensate. If this happens, populations at states with lower winding numbers, namely $`m=1`$ and $`m=0`$ states, start to increase. This implies that angular momentum will decrease. Physically this is nonsense since it is clear that angular momentum per particle must increase monotonically with the rotation frequency. Consequently we expect an upper critical frequency $`\stackrel{~}{\mathrm{\Omega }}_{up}`$ when the coreless vortex ceases to be the ground state of the system.

We can get a rough estimate for this upper limit by assuming that the slope depends linearly on $`\stackrel{~}{\mathrm{\Omega }}`$ and extrapolating the behavior to higher rotation frequencies. For the parameters used before this gives us an upper critical frequency $`\stackrel{~}{\mathrm{\Omega }}_{up}0.15`$, which is again in good agreement with the solutions of the GP equations .

The ground state at higher rotation frequency is more complicated. The superfluid velocity field slightly above the upper critical rotation frequency is still rotationally symmetric, but the order parameter is no longer a simple mixture of components with winding numbers $`0`$, $`1`$, and $`2`$. In Figure 4 we show an example of the order parameter when $`\stackrel{~}{\mathrm{\Omega }}=0.17`$. In the $`m=0`$ component we can see a regular array of $`4`$ vortices. The $`m=\pm 1`$ components fill the vortex cores. Although it is not clear from the figure, the $`m=\pm 1`$ components have a nonzero angular momenta. The velocity field of this configuration is rotationally symmetric and qualitatively resembles that in Figure 1.

In reality $`\beta (r)`$ is not linear, but its slope vanishes close to the edge of the condensate. It seems that the maximum boundary value for the bending angle is $`\beta _{max}=\frac{3}{4}\pi `$ and AT-vortex is never the ground state of the system. We believe that this, somewhat surprising, boundary value is caused by a subtle effect due to kinetic energy. Ordinarily kinetic energy of the BEC can be safely ignored, but close to the boundary of the condensate where gradients are large its effects are pronounced. We conjecture that the coreless vortex structure with $`\beta _{max}>\frac{3}{4}\pi `$ is dynamically unstable and that this instability is due to the physics close to the condensate edge. But we freely admit that at the time of writing this paper, a detailed understanding of the spinor condensate dynamics close to the boundary of the condensate is still lacking.

## IV Antiferromagnetic ground state

It should be understood that the results change dramatically with the sign of $`\lambda _a`$ . While a coreless vortex can carry an infinitesimal amount of angular momentum, this state is not energetically favored if $`\lambda _a>0`$. For polar condensates the ordinary non-rotating ground state $`\zeta =(0,1,0)^T`$ remains the ground state until a certain critical frequency is reached. This critical frequency is somewhat higher than the corresponding critical frequency for a (ferromagnetic) coreless vortex, but well below the critical frequency for the stabilization of an ordinary singly quantized vortex.

Just above this threshold the ground state of the system is a state with angular momentum $`1/2`$ per particle (see Figure 5). This state is locally magnetized and cannot be represented as a local rotation of the non-rotating ground state. This is easily understood, since there is no such local rotation of the polar state that results in a texture having non-zero angular momenta. While the state is locally magnetized, the number of atoms at $`m=\pm 1`$ states is still the same and the system is globally non-magnetized. We can choose a convenient basis by a global rotation of the numerically calculated wavefunction. Such a procedure reveals that this state can be intepreted as a state where population is evenly split between $`m=\pm 1`$ states. Another one of these states carries a vortex and the other one fills the vortex core. Incidentally, the dynamical stability of this structure was predicted in Ref. . Here we have shown that within a certain range of rotation frequencies, this state is also the minimum energy configuration. Very interesting point is that the winding number $`1/2`$ vortex of the polar spinor condensate can viewed as the BEC counterpart to the Alice string familiar from the particle physics .

Even though the behavior of the system depends strongly on the sign of $`\lambda _a`$, it does not depend strongly on its magnitude. $`\lambda _a`$ could be increased threefold without large qualitative or quantitative changes. Also, the number of particles could be changed by, at least, a factor of two without important qualitative changes. The range of rotation frequencies when this state is the ground state is quite narrow. In our simulations this range was only on the order of $`0.01`$. If the rotation frequency is beoynd this narrow range, the ground state is a state where population is evenly split between $`m=\pm 1`$ states, but in this new state both components carry a vortex. These vortices are not overlapping and the total angular momentum per particle is less than $`1`$. These results are analogous to the results for scalar condensates, where it has been shown that in a parabolic trap a multiply quantized vortex is never the ground state . In a parabolic trap it is always energetically favorable to break a multiply quantized vortex into pairs of singly quantized vortices.

## V Magnetic field

So far we have assumed that the chemical potentials of the different $`m`$-states are equal. In this section we discuss the effect of an external magnetic field pointing along the rotation axis. If the magnetic field has components also in $`x`$\- and $`y`$-directions, things can be different, as will be discussed in the next section.

External magnetic field will shift the chemical potential of the component $`m`$ into $`\mu _m=\mu _0+mE_Z`$, where $`E_Z`$ is the Zeeman shift and $`\mu _0`$ is the chemical potential of the $`m=0`$ component. An external magnetic field will tend to magnetize the system and move more population into the state with the lowest chemical potential. It is simple to show, that for a non-rotating homogeneous spinor condensate the Zeeman shift required for a complete magnetization in the direction of the external magnetic field is on the order of $`\lambda _an`$, where $`n`$ is the density . Note that the magnetizing field is very small e.g. for sodium atoms with reasonable densities it is about $`B=10^410^3`$ G.

Assume a slowly rotating system prepared in the absence of the magnetic field and allowed to equilibriate into the ground state. We then turn on the magnetic field which is aligned with respect to the rotation axis and the angular momentum of the system. With an increasing magnetic field the ground state of the system evolves towards the ground state of the scalar condensate in $`m=1`$ or $`m=1`$ component. Which component is favored depends on the direction of the magnetic field. As the rotation frequency is small, this ground state has a vanishing angular momentum. Therefore angular momentum of the ground state decreases with increasing magnetic field.

The reverse process of starting with a large magnetic field along the $`z`$-direction and then turning it off behaves differently. The ground state of the scalar GP equation is metastable under such change. Initially with a large magnetic field the ground state is essentially a scalar condensate without angular momentum. It is hard to relax such a state into a coreless vortex state with finite angular momentum. This behavior is reminescent of the hysteresis for the vortex nucleation in scalar condensates.

## VI Creation

As it was already mentioned, the behaviour of the system may be different if the external magnetic field has other non-zero components than the one pointing in $`z`$-direction. This is dramatically shown in the topological vortex formation  which was recently demonstrated experimentally . In this experiment a magnetic field configuration of the Ioffe-Pritchard trap was changed in time. In particular the $`z`$-component of the magnetic field was inverted. Once the $`z`$-component is inverted, the condensate wavefunction has acquired a winding number $`2`$. It should be noted that this process continuously transforms a condensate with vanishing angular momenta into a condensate with an angular momentum $`2`$ per particle. When the spreading of the wavefunctions can be ignored this process is equivalent with coreless vortex transformations. In this respect coreless vortices have already been experimentally created. In a Ioffe-Pritchard magnetic trap a coreless vortex is not stable since all hyperfine states are not trapped. The stabilization of the coreless vortex structures would require an optical trap.

## VII Conclusions

We have studied the ground state properties of rotating spin-1 condensates. We have found out that for a sufficiently small rotating frequency the ground state of a ferromagnetic spinor condensate is a coreless vortex. However, because of the finite size of a condensate a certain critical rotation frequency has to be exceeded for a coreless vortex to appear. We have also investigated antiferromagnetic condensates which in contrast to coreless vortices have rotating ground states exhibiting the discrete staircase behavior.

All the simulations presented in this paper have been performed in two dimensions. Consequently they are expected to describe well pancake shaped condensates. We have also performed simulations with full three-dimensional Gross-Pitaevskii equations and convinced ourselves that our results are not specific to two dimensions. While a detailed study in three dimensions is at present too time-consuming, it does seem that our results apply qualitatively even in a cigar shaped trap.

In this paper we have focused solely on the ground state. Naturally the ground state is one thing, but getting the system into it in a reasonable time is quite another. The energy landscape of the rotating spinor condensate is complex and metastable states abound. This may cause problems in the experiments as the relaxation time to the true ground state can become very long. Also the problem of nucleation of coreless vortices may cause problems, especially if the system is initially prepared as a pure stationary condensate. Nevertheless with some effort all these complications can be overcome, and the ground state of the rotating spinor condensate promises to be an exciting field of research.

## VIII Acknowledgements

Authors acknowledge the Academy of Finland (project 50314) for financial support. In addition J.-P. M. was supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is supported by the Nederlandse Organisatie voor Wetenschaplijk Onderzoek (NWO). A. C. thanks the Magnus Ehrnrooth foundation for financial support.