text:
# 1 Introduction

## 1 Introduction

In a range of applications we often seek curves and surfaces that have an aesthetically pleasing ‘roundedness’ to them. One way to make this concept concrete is through looking for globally-optimised ‘minimum variation curves’ . The philosophy behind this idea is straightforward. We usually prefer curves that are close to circular over curves with sharp turns. This is particularly true when designing camera trajectories, where sudden changes in curvature can have a very disorienting effect. Circular paths are characterised by having constant curvature, so a natural idea in forming interpolations between control points is to find a curve that minimises the total change in curvature. The problem with such a strategy is that these curves can be extremely hard to compute. If one adopted a variational strategy, with endpoint conditions, the equations for the curve can be as high as fifth order and are even more difficult to treat than those of elasticity. Such equations can only be solved numerically and do not have straightforward, controllable, analytic solutions. The problem is even more acute if multiple control points are involved, as even numerical computation can be extremely difficult.

A more straightforward, local scheme that provides smooth interpolations was introduced by Wenz  and later extended by Szilvási-Nagy & Vendel . The idea explored by these authors is to generate curves that are as close as possible to circles. Given four points $`X_0,\mathrm{},X_3`$ we construct the circles $`C_1`$ through $`X_0`$, $`X_1`$ and $`X_2`$, and $`C_2`$ through $`X_1`$, $`X_2`$ and $`X_3`$. The curve between $`X_1`$ and $`X_2`$ is then formed by smoothly interpolating between the two circles. This idea was further extended by Séquin & Yen  and Séquin & Lee , who introduced an angle-based circle blending scheme. The angle-based scheme gives better results than the earlier, midpoint scheme, and we argue here that it is geometrically the ‘correct’ one. Séquin & Lee also showed how to achieve $`G_2`$-continuity (and higher order continuity, if desired), and demonstrated the value of angle-based blending for interpolating over the surface of a sphere.

In this paper we further explore the geometry associated with circle blending, following the methods developed by Séquin and his coworkers. The essential idea is that the natural way to transform between two circles is via a conformal transformation. Conformal transformations leave angles invariant, but can alter distances. Euclidean transformations are the subset of conformal transformations that also leave distances invariant. Conformal transformations in a plane can take any three chosen points to any three image points. As such, they can transform a line or circle into any other circle. In this geometry, straight lines are examples of circles that pass through the point at infinity. By exploiting the features of conformal geometry, we can write robust code that treats (straight) lines and circles in a single, unified manner. This eliminates the need to check for special cases. Similarly, in three dimensions, planes and spheres are treated as examples of the same object. So a single routine can interpolate between points on a sphere, and will reduce to the planar case when four points happen to lie in a plane.

To fully exploit the advantages of conformal geometry we work in the mathematical framework of geometric algebra . This algebra treats points, lines, circles, planes and spheres, and the transformations acting on them, in a unified algebraic framework. A number of authors have argued for the advantages of the conformal geometric algebra framework for computer graphics applications . The present work should be viewed in this context. We show how complex problems such as finding the conformal transformation between a line and a circle reduce to simple, robust expressions in the geometric algebra framework. As a further application we show how the same framework naturally extends to sphere-blending over a surface. This suggests a new method of characterising surfaces that does not require the concept of swept curves.

This paper starts with an introduction to conformal geometric algebra. This introduction is self-contained, but to keep its length down a number of concepts are introduced with a minimum of explanation. We then turn to the question of how to mathematically encode transformations between circles. We find the conformal transformation that achieves this and explore its properties. Some subtleties involving the orientation of the transformation are explained, and we demonstrate how they are easily resolved in the conformal framework. We then provide a series of examples of blended curves, and illustrate the effects of demanding higher-order $`G`$-continuity. We finish by introducing a method of sphere blending and discuss the potential of this idea for encoding surfaces.

## 2 Conformal geometry and Euclidean space

The starting point for our description of geometry is the conformal group. This marks a radical departure from conventional descriptions of Euclidean space based on projective geometry and homogeneous coordinates. The main advantage of basing the description in a conformal setting is that distance is encoded simply, making the geometry well suited to describing the real three-dimensional world.

Suppose we start with an $`n`$-dimensional Euclidean vector space $`^n`$. The conformal group consists of the set of all transformations of $`^n`$ that leave angles invariant. These include translations and rotations, so the conformal group includes the set of Euclidean transformations as a subgroup. The conformal group on $`^n`$ has a natural representation in terms of rotations in a space two dimensions higher, with signature $`(n+1,1)`$. So, in the same way that projective transformations are linearised by working in a space one dimension higher than the Euclidean base space, conformal transformations are linearised in a space two dimensions higher. The conformal representation of points in Euclidean three-space consists of vectors in a 5-dimensional space. While this may appear to be an unnecessary abstraction, working in this five-dimensional space does bring a number of advantages.

To exploit the conformal representation we need a standard representation for a Euclidean point in the five-dimensional conformal space. Given that we are occupying a space two dimensions higher, two constraints are required to specify a unique point. The first of these is that our underlying representation is homogeneous, so $`X`$ and $`\lambda X`$ represent the same point in Euclidean space. The second constraint is that the vector $`X`$ is null,

$$X^2=0.$$
(1)
(The existence of null vectors is guaranteed by the fact that the conformal vector space has mixed signature.) This is essentially the only further constraint that can be enforced which is consistent with homogeneity and invariant under orthogonal transformations in conformal space. Now suppose that $`e_1`$, $`e_2`$ and $`e_3`$ represent three vectors in the three-dimensional base space, and we add to these the vectors $`e_0`$ and $`e_4`$. These satisfy

$$e_0^2=1,e_4^2=+1,$$
(2)
and all 5 vectors $`\{e_0,\mathrm{},e_4\}`$ are orthogonal. From the two extra vectors we define the two null vectors $`n`$ and $`\overline{n}`$ by

$$n=e_4+e_0,\overline{n}=e_4e_0,n\overline{n}=2.$$
(3)
(It is a straightforward exercise to confirm that these two vectors both have zero magnitude.) From these we need to chose a vector to represent the origin. This is conventionally taken as $`\frac{1}{2}\overline{n}`$. The vector $`X`$ can therefore be written as

$$X=2x\overline{n}+\alpha n$$
(4)
where $`x`$ is the equivalent three-vector representation of the point in Euclidean space. The variable $`\alpha `$ must be chosen so that $`X`$ is null, which fixes $`\alpha =x^2`$. We therefore arrive at the following representation of a point in conformal space:

$$X=2x+x^2n\overline{n}.$$
(5)
This representation can be arrived at using more geometrical reasoning by considering a stereographic projection . It is immediately clear from equation (5) that $`n`$ represents the point at infinity.

The Euclidean coordinates of the point $`x`$ are recovered from $`X`$ via the homogeneous relation

$$x_i=\frac{Xe_i}{Xn},i=1,2,3.$$
(6)
This relation confirms that $`X`$ and $`\lambda X`$ define the same Euclidean point. In principle, one could allow for $`\lambda `$ to be negative, but in practice this should be avoided. Restricting to positive $`\lambda `$ is is equivalent to stating that $`X`$ must always have a positive $`e_0`$ component. The reason for this restriction can be seen in figure 1. The equation $`X^2=0`$ generates a cone structure, and we only want to employ one half of this to represent points. Maintaining a positive $`e_0`$ component throughout all algorithms enables us to keep track of orientations consistently.

Given a point $`a`$, represented by the unnormalised vector $`A`$, we can place $`A`$ in the standard form of equation (5) with the map

$$A2\frac{A}{An}.$$
(7)
One should be wary of employing this map when writing code, as the right-hand side is singular if $`A`$ happens to be the point at infinity. As is the case for projective geometry, it is better practice to let the normalisation run free, and only use equation (6) in the final stage to recover the coordinates.

The power of the conformal representation starts to become clear when we consider the inner product of points. Suppose that $`X`$ and $`Y`$ are the conformal representations of the points $`x`$ and $`y`$ respectively, both in the standard form of equation (5). Their inner product is

$`XY`$ $`=\left(x^2n+2x\overline{n}\right)\left(y^2n+2y\overline{n}\right)`$
$`=2x^22y^2+4xy`$
$`=2(xy)^2.`$ (8)
This is the essential result that underpins the conformal approach to Euclidean geometry. The inner product in conformal space encodes the distance between points in Euclidean space. This is the reason why points are represented with null vectors — the distance between a point and itself is zero.

Generalising equation (8) to unnormalised vectors, the distance between points can be written

$$|xy|^2=2\frac{XY}{XnYn}.$$
(9)
Any transformation of $`X`$ and $`Y`$ that leaves this product invariant must represent a symmetry of Euclidean space. Because such a transformation must leave the inner product invariant, it must be an orthogonal transformation in conformal space. Such a transformation ensures that null vectors remain null vectors, and so continue to represent points. The transformation must also leave $`n`$ invariant if the distance is to remain unchanged. We can now see that the Euclidean group is the subgroup of conformal transformations that leaves invariant the point at infinity. In what follows we will find applications for both Euclidean and more general conformal transformations.

## 3 Conformal Geometric Algebra

Geometric algebra was introduced by the nineteenth century mathematician W.K. Clifford, and has found many applications in physics . Often in this work the name Clifford algebra is used, but for applications in geometry it is becoming increasingly common to see Clifford’s original name of geometric algebra to describe the field. A geometric algebra is constructed over a vector space with a given inner product. The geometric product of two vectors $`a`$ and $`b`$ is defined to be associative and distributive over addition, with the additional rule that the square of any vector is a scalar,

$$aa=a^2.$$
(10)
If we write

$$ab+ba=(a+b)^2a^2b^2$$
(11)
we see that the symmetric part of the geometric product of any two vectors is also a scalar. This defines the inner product, and we write

$$ab=\frac{1}{2}(ab+ba).$$
(12)
The remaining, antisymmetric part of the geometric product returns the outer or exterior product familiar from projective geometry. We write this as

$$ab=\frac{1}{2}(abba).$$
(13)
The totally antisymmetrised sum of geometric products of vectors defines the exterior product in the algebra.

The geometric product of two vectors can now be written

$$ab=ab+ab.$$
(14)
Under the geometric product, orthogonal vectors anticommute and parallel vectors commute. The product therefore encodes the basic geometric relationships between vectors. Now that we know how to multiply vectors together, it is straightforward to construct the entire geometric algebra of a given vector space. This is facilitated by introducing an orthonormal frame of vectors $`\{e_i\}`$. These satisfy

$$e_ie_j=\eta _{ij},$$
(15)
where $`\eta _{ij}`$ is the metric tensor. For a space of signature $`(p,q)`$, $`\eta _{ij}`$ is a diagonal matrix consisting of $`p`$ $`+1`$s and $`q`$ $`1`$s. The space of interest to us here is the conformal space of signature $`(4,1)`$. A basis for this is provided by the vectors $`e_0,\mathrm{},e_4`$, with $`e_0`$ having negative square. The algebra generated by these vectors has 32 terms in total, and is spanned by

$$\begin{array}{ccccccc}& 1& \{e_i\}& \{e_ie_j\}& \{e_ie_je_k\}& \{Ie_i\}& I\\ \hfill \text{grade}& 0& 1& 2& 3& 4& 5\\ \hfill \text{dimension}& 1& 5& 10& 10& 5& 1.\end{array}$$

We refer to this algebra as $`𝒢(4,1)`$. The term ‘grade’ is used to refer to the number of vectors in any exterior product. The dimensions of each graded subspace are given by the binomial coefficients.

The highest grade term in $`𝒢(4,1)`$ is called the pseudoscalar and is given the symbol $`I`$. This is defined by

$$I=e_0e_1e_2e_3e_4.$$
(16)
The pseudoscalar commutes with all elements in the algebra, a feature of odd-dimensional algebras, and the $`(4,1)`$ signature of the space implies that the pseudoscalar satisfies

$$I^2=1.$$
(17)
So, algebraically, $`I`$ has the properties of a unit imaginary. But it also plays a definite geometric role, as multiplication by the pseudoscalar performs a duality transformation in conformal space. A matrix representation of $`𝒢(4,1)`$ can be constructed in terms of $`4\times 4`$ complex matrices. These can be found in the physics literature in the guise of the Dirac matrices . For practical applications this matrix representation has little value and one is better off coding up the algebraic rules explicitly.

A general element of $`𝒢(4,1)`$ is called a multivector and can consist of a sum of terms all grades in the algebra. Arbitrary elements of $`𝒢(4,1)`$ can be added and multiplied together. A multivector that consists only of terms of a single grade is said to be homogeneous. The geometric product of a pair of homogeneous multivectors decomposes as follows:

$$A_rB_s=A_rB_s_{|rs|}+A_rB_s_{|rs|+2}+\mathrm{}+A_rB_s_{r+s}.$$
(18)
The angle brackets $`M_r`$ are used to denote the projection onto the grade-$`r`$ terms in $`M`$. The dot and wedge symbols are used to generalise the inner and outer products to the lowest and highest grade terms in a general product:

$`A_rB_s`$ $`=A_rB_s_{|rs|}`$
$`A_rB_s`$ $`=A_rB_s_{r+s},`$ (19)
where $`A_r`$ and $`B_s`$ are homogeneous multivectors of grade $`r`$ and $`s`$ respectively. One only ever uses the dot and wedge symbols when applied to homogeneous multivectors. Of particular importance are the inner and outer products with a vector. These satisfy

$`aA_r`$ $`=\frac{1}{2}(aA_r(1)^rA_ra)`$
$`aA_r`$ $`=\frac{1}{2}(aA_r+(1)^rA_ra).`$ (20)
Some straightforward algebra confirms that the outer product is associative .

An important algebraic operation applied to multivectors is reversion. This plays an analogous role to transposition in matrix algebra. The reverse of a multivector $`M`$ is denoted by $`\stackrel{~}{M}`$ and is defined by reversing the order of all geometric products of vectors in $`M`$. An arbitrary multivector in $`𝒢(4,1)`$ can be written as

$$M=\alpha +a+B+T+Ib+I\beta $$
(21)
where $`\alpha `$ and $`\beta `$ are scalars, $`a`$ and $`b`$ are vectors, $`B`$ is a bivector and $`T`$ is a trivector. The reverse of $`M`$ is then given by

$$\stackrel{~}{M}=\alpha +aBT+Ib+I\beta .$$
(22)
A conformal transformation in the Euclidean space $`^3`$ is represented by an orthogonal transformation in $`𝒢(4,1)`$. Of particular relevance here are special conformal transformations. These have determinant $`+1`$, and correspond to orientation-preserving transformations on $`^3`$. In geometric algebra a special conformal transformation, applied to an arbitrary multivector $`A`$, can be written in the compact form

$$ARA\stackrel{~}{R}.$$
(23)
Here $`R`$ is an even-grade element (grades 0, 2 and 4) satisfying

$$R\stackrel{~}{R}=1.$$
(24)
Because our representation is homogeneous, this normalisation constraint can usually be relaxed to allow for $`R\stackrel{~}{R}`$ to be an arbitrary positive scalar. The element $`R`$ is called a rotor. As an example, the (unnormalised) rotor for the Euclidean translation between $`x`$ and $`y`$ can be written

$$T_{xy}=(nY)nX+(nX)Yn,$$
(25)
where $`X`$ and $`Y`$ are the conformal equivalents of $`x`$ and $`y`$ respectively. Every rotor can be written in term of the exponential of a bivector as

$$R=\pm \mathrm{exp}(B/2).$$
(26)
The bivector $`B`$ is the generator of the rotation, and is an element of the Lie algebra of the conformal group.

## 4 Geometric primitives in conformal algebra

We have now established how points in $`^3`$ are represented by null vectors in the conformal algebra $`𝒢(4,1)`$. We now need to establish how higher dimensional objects are represented. Given a null vector $`Y`$, one view of the point this represents is as the solution space of the equation

$$YX=0,X^2=0.$$
(27)
where here $`X`$ is treated as a vector variable. This is a similar idea to projective geometry, and generalises immediately to higher-grade objects. Given a multivector $`A_r`$, the outer product of $`r`$ distinct vectors, the geometric object that this represents is the solution space of

$$A_rX=0,X^2=0.$$
(28)
This pair of equations is clearly homogeneous in $`X`$. If a conformal transformation is applied to $`A_r`$, taking it to $`RA_r\stackrel{~}{R}`$, then the solution space transforms in exactly the same way (that is, $`XRX\stackrel{~}{R}`$). In this manner conformal transformations are easily applied to higher-grade objects.

After vectors, the next objects to consider are bivectors. Given a bivector $`B`$ (formed from the outer product of two vectors), the solution space of equation (28) depends on the sign of $`B^2`$. If $`B^2`$ is negative there are no solutions, if $`B^2=0`$ there is one solution, and if $`B^2`$ is positive there are two solutions. This final case corresponds to the situation where $`B=XY`$, where $`X`$ and $`Y`$ are a pair of null vectors. Given a bivector of this form it is straightforward to recover $`X`$ and $`Y`$ , so a bivector can encode a pair of points.

Next consider a trivector $`L`$. The main example of relevance is the trivector formed from the outer product of three points $`X_1`$, $`X_2`$ and $`X_3`$, so

$$L=X_1X_2X_3.$$
(29)
Any vector $`X`$ that is a solution of $`LX=0`$ must be a linear combination of $`X_1`$, $`X_2`$ and $`X_3`$. The fact that $`X`$ is null and has arbitrary scale implies that in $`^3`$ the solution space is one-dimensional. To see what this space is, consider the null vector $`C`$, where

$$C=LnL.$$
(30)
If $`X`$ satisfies $`XL=0`$, then $`X`$ commutes with $`L`$. It follows that

$$\frac{XC}{XnCn}=\frac{L^2}{Cn}=\text{constant}.$$
(31)
The solution set $`X`$ therefore consists of points at an equal distance from some point $`C`$. It follows that the trivector $`L`$ represents the circle through $`X_1`$, $`X_2`$ and $`X_3`$. If this circle also passes through the point at infinity, the line is straight. That is, the test that a line $`L`$ is straight is that

$$Ln=0.$$
(32)
Similarly, the straight line through $`X_1`$ and $`X_2`$ is given by

$$L=X_1X_2n.$$
(33)
Conformal algebra therefore treats lines and circles in a unified manner as trivectors in $`𝒢(4,1)`$. This is to be expected, as a conformal transformation can always transform a line into a circle.

The radius $`\rho `$ of the circle $`L`$ is found from

$$\rho ^2=\frac{L^2}{(Ln)^2}.$$
(34)
So again we see how a metric object (the radius) is recovered from a homogeneous representation of a geometric entity (the circle $`L`$). The angle $`\theta `$ between two circles (or lines) is found from

$$\mathrm{cos}\theta =\frac{L_1L_2}{|L1||L2|}.$$
(35)
As expected, this is invariant under the full conformal group.

Similar considerations apply to the 4-vector $`S`$ defined by the four points $`X_1,\mathrm{},X_4`$,

$$S=X_1X_2X_3X_4.$$
(36)
The object defined by $`S`$ is the unique sphere through the four points (which cannot all lie on a line or circle if $`S0`$). The centre of the sphere is given by $`SnS`$, and if the centre lies at infinity the sphere reduces to a plane. As with the case of a straight line, the plane defined by the three points $`X_1`$, $`X_2`$ and $`X_3`$ is

$$P=X_1X_2X_3n.$$
(37)
Lines, circles, planes and spheres can all be intersected in a straightforward manner in conformal geometric algebra. The two cases of greatest significance are those of a circle and a sphere, and of two spheres. Treating the latter first, any two spheres intersect in a circle (provided they touch). This reduces to a line if both spheres are in fact planar. The intersection $`L`$ is found simply by

$$L=(IS_1)S_2,$$
(38)
which is a trivector, as required. Similarly, the intersection of a circle $`C`$ and a sphere $`S`$ is given by

$$B=(IC)S.$$
(39)
This is a bivector, as a circle can intersect a sphere or plane in zero, one, or two points. In the case where $`C`$ is a straight line and $`S`$ is a flat plane, one of the points of intersection is at infinity. As we proceed we will see how the orientations implicit in lines and surfaces are neatly encoded in these intersection formulae.

## 5 Circle blending

Suppose we are given a series of $`n+1`$ points $`x_0,\mathrm{},x_n`$ and we seek a smooth curve through all of these points. The idea behind the circle blending scheme is to find a curve that is as close as possible to a circle and which passes through all of the points. A typical setup is illustrated in figure 2. Consider the four points $`x_{i1}\mathrm{}x_{i+2}`$. We form the circle $`C_i`$ through $`x_{i1}`$, $`x_i`$ and $`x_{i+1}`$, and the circle $`C_{i+1}`$ through $`x_i`$, $`x_{i+1}`$ and $`x_{i+2}`$. In the region between $`x_i`$ and $`x_{i+1}`$ we need a curve that blends smoothly between the circles $`C_i`$ and $`C_{i+1}`$.

A range of schemes exist for interpolating between circles, but one is naturally picked out from the viewpoint of conformal geometry. We are presented with two circles $`C_1`$ and $`C_2`$, both of which share two points in common:

$`C_1`$ $`=X_0X_1X_2`$
$`C_2`$ $`=X_1X_2X_3.`$ (40)
Given the (conformal) points $`X_1`$, $`X_2`$, and the circles $`C_1`$ and $`C_2`$, there are no further objects to specify. First we need to define family of circles between $`C_1`$ and $`C_2`$. This is easily achieved by finding the transformation between them. The transformation must map between two circles, so is a conformal transformation. The only bivector generator for such a transformation is the antisymmetrised product between $`C_1`$ and $`C_2`$. A transformation governed by this generator is simply a rotation in $`𝒢(4,1)`$. If we normalise the circles by defining

$$\widehat{C_i}=\frac{C_i}{|C_i|}$$
(41)
then the interpolated circle is simply

$$\widehat{C}_{12}(\lambda )=\frac{1}{\mathrm{sin}(\theta )}\left(\mathrm{sin}\left((1\lambda )\theta \right)\widehat{C}_1+\mathrm{sin}(\lambda \theta )\widehat{C}_2\right).$$
(42)
Here $`\theta `$ is the angle between the circles, defined by

$$\mathrm{cos}(\theta )=\widehat{C}_1\widehat{C}_2.$$
(43)
(The case of $`\theta >\pi `$ will be addressed shortly.) This method of interpolating between circles recovers the angle-blending scheme of Séquin & Yen , thus providing a firm geometric reason for preferring that scheme. At this point there is no reason for $`C_1`$ and $`C_2`$ to lie in the same plane. If they lie on different planes then the four points $`X_0,\mathrm{},X_3`$ define a sphere. Since the interpolated circles can be written

$$\widehat{C}_{12}(\lambda )=\frac{1}{\mathrm{sin}(\theta )}X_1X_2\left(\mathrm{sin}\left((1\lambda )\theta \right)\frac{X_0}{|C_1|}+\mathrm{sin}(\lambda \theta )\frac{X_2}{|C_2|}\right),$$
(44)
we see that all points on these circles are combinations of $`X_0\mathrm{}X_3`$. It follows that all of the intermediate circles also lie of the sphere defined by $`X_0,\mathrm{},X_3`$, which is a desirable property for a range of applications . The basic interpolation scheme is illustrated in figure 3

Now that we have a straightforward means of encoding the circle blends, we need to parameterise the actual trajectory between the points. A simple means for achieving this is to start with the straight line between $`X_1`$ and $`X_2`$. This path is parameterised by

$$Y_{12}(\lambda )=(1\lambda )X_2nX_1\lambda X_1nX_2+\lambda (1\lambda )X_1X_2n,$$
(45)
and the line itself is described by

$$L_{12}=X_1X_2n.$$
(46)
The rotor that transforms between the line $`L_{12}`$ and the circle $`C_{12}(\lambda )`$ is simply

$$R_{12}(\lambda )=1+\widehat{C}_{12}(\lambda )\widehat{L}_{12}.$$
(47)
So once the (normalised) blended circle $`\widehat{C}_{12}(\lambda )`$ is found, the path itself is given by transforming from the straight line to the circle. The path between $`X_1`$ and $`X_2`$ is therefore simply

$$X_{12}(\lambda )=R_{12}(\lambda )Y_{12}(\lambda )\stackrel{~}{R}_{12}(\lambda ).$$
(48)
This is extremely simple to code up. We never have to calculate the centre of a circle, and the algorithm deals equally easily with straight lines or circles. The algorithm keeps all intermediate points on the sphere defined by each group of four points. If all of the base points lie on the same sphere, then we compute smooth curves over this sphere. An example of this is shown in figure 4.

## 6 Problems and special cases

Provided the sequence of any four points are unique, the circles to be blended and the intermediate straight line are all well defined. The only aspects of the algorithm that can be problematic are concerned with the definition of the angle $`\theta `$ in equation (42). There are two cases over which care must be taken, both of which can be treated fairly easily. The first obvious problem is that the definition of the blended segment breaks down if $`\mathrm{sin}(\theta )=0`$, which occurs when $`\mathrm{cos}(\theta )=+1`$ and when $`\mathrm{cos}(\theta )=1`$. The former case corresponds to $`\theta =0`$ and implies that the two curves are the same. This is easily caught by setting a threshold value for $`\theta `$ below which eqn (42) is replaced by the small-$`\theta `$ limit of

$$C_{12}(\lambda )=(1\lambda )\widehat{C}_1+\lambda \widehat{C}_2.$$
(49)
The case of $`\theta =\pi `$ is more complicated, and corresponds to a somewhat pathological example. This occurs when the two circles only differ in their orientation, so $`\widehat{C}_2=\widehat{C}_1`$. A sample configuration is shown in figure 5. All four points lie in a plane, so the interpolated curve must also lie in this plane. It is possible to define an interpolation scheme for this case. If we let $`P`$ denote the plane containing the circles, then the generator of the transformation must act in the $`P`$ plane, leaving $`X_1`$ and $`X_2`$ invariant. The bivector generator for the transformation is therefore

$$B_{12}=(X_1X_2)P.$$
(50)
This produces the type of interpolated curve shown in figure 5. The problem here is that this case is not smoothly connected to the general case. Imagine moving one of the control points in figure 5 slightly out of the page. The four points then define a sphere, and the interpolated curve will lie on this sphere. This is a totally different curve to the planar case. If we are only interested in planar plots, then everything is well-behaved and the above case can be easily dealt with. But if we are interested in curves in three-dimensional space, this case is best avoided altogether with the addition of further control points.

The case of oppositely oriented circles relates to the second key issue, which is what happens if the circles are greater than $`180^{}`$ apart? The arccos function will not return the correct angle, and the interpolation scheme will effectively blend the wrong way. This problem has an elegant solution within the conformal framework. The question reduces to finding the correct midpoint circle $`C_{1,2}`$. This is given by

$$C_{1,2}=\{\begin{array}{cc}\widehat{C}_1+\widehat{C}_2,\hfill & \theta <\pi \hfill \\ (\widehat{C}_1+\widehat{C}_2),\hfill & \theta >\pi \hfill \end{array}.$$
(51)
The case of $`\theta =\pi `$ corresponds to equal and opposite circles, which is the one awkward case we have to avoid. All we need is a test to decide which is the correct mid-circle for the geometry, and the arccos function can be used to determine $`\theta /2`$ unambiguously from

$$\mathrm{cos}(\theta /2)=\widehat{C}_1\widehat{C}_{1,2}.$$
(52)
To find a suitable test we consider the plane of points equidistant from $`X_1`$ and $`X_2`$. This plane is defined by

$$\mathrm{\Pi }=I(X_1X_2)n.$$
(53)
All planes are oriented, so the sign here is important. This choice corresponds to the normal to the plane (defined by the right-hand rule) pointing in the direction from $`X_1`$ to $`X_2`$. Now suppose that we intersect a circle $`C`$ with the plane $`\mathrm{\Pi }`$. The result is two points encoded in the bivector $`B`$, where

$$B=I(C\mathrm{\Pi }+\mathrm{\Pi }C).$$
(54)
The bivector $`B`$ contains the points in the order

$$B=X_fX_i,$$
(55)
where $`X_f`$ is the point where the circle intersects the plane from the negative to positive side, and $`X_i`$ is the opposite point (see figure 6).

Now suppose that $`M_1`$ and $`M_2`$ are the midpoints of the circle segments $`C_1`$ and $`C_2`$ respectively. Both lie in the plane $`\mathrm{\Pi }`$. The correct mid-circle should have its $`X_f`$ intersection point closer to the midpoint of $`M_1`$ and $`M_2`$ than its $`X_i`$ point. If this is not the case, then the circle has the wrong orientation. We therefore only need form the quantity

$$\alpha =\left((M_1nM_2+M_2nM_1)n\right)B$$
(56)
where $`B`$ is given by equation (54) with $`C`$ defined by

$$C=\widehat{C}_1+\widehat{C}_2.$$
(57)
The test is now

$$C_{1,2}=\{\begin{array}{cc}\widehat{C}_1+\widehat{C}_2,\hfill & \alpha <0\hfill \\ (\widehat{C}_1+\widehat{C}_2),\hfill & \alpha >0\hfill \end{array}.$$
(58)
This conditional statement covers all special cases (where various circles are degenerate or reduce to straight lines) and the border cases of $`\alpha =0`$ is the single degenerate case mentioned above. With the mid-circle found, the angle $`\theta `$ is computed, and the interpolation scheme of equation (42) can be followed.

A further problem with the blending scheme is highlighted in figure 7. The geometry of the control points ensures that an unwanted inversion arises in the middle of the plot. This is because the gradient at each control point is fixed entirely by the two points on either side (see section 7), and so in this case is forced to be flat, as illustrated. The problem can be removed in a number of ways. One of the simplest is to define a series of midpoints, and then repeat the blending algorithm. The result of a single iteration of this type is also shown in figure 7. Of course, one can go further and continue to recursively define midpoints, which will also generate a smooth blend. This latter approach has the attractive feature that it removes the need for any trigonometric evaluations, as the mid-circle is found easily using the method described above.

## 7 Continuity

An important issue in any interpolation scheme is the order of continuity at the control points. In this section we will only consider $`G`$-continuity. To achieve $`C`$-continuity some reparameterisation will have to occur on each spline. For the case of a camera fly-by, this can be achieved by specifying the desired velocity along the curve.

Our blending scheme is defined by equation (48), and before we can differentiate this we need a result for the derivative of a rotor. The definition of a (normalised) rotor is sufficient to prove that we can always write

$$\frac{R}{\lambda }=\frac{1}{2}RB$$
(59)
where $`B`$ is a bivector, which in general will also be a function of $`\lambda `$. (The factor of $`1/2`$ is a useful convention). To find the tangent vector to the (conformal) curve $`X(\lambda )`$ we form the line

$$T=XX^{}n$$
(60)
where the dash denotes the derivative with respect to $`\lambda `$. To see how this behaves we return to equation (48) and form the derivative at $`\lambda =0`$. For convenience we will assume that the rotor in equation (47) has been normalised, though this is unimportant when forming the tangent vector. In evaluating the derivatives we can use the results that

$$R(\lambda )X_1\stackrel{~}{R}(\lambda )=X_1,R(\lambda )X_2\stackrel{~}{R}(\lambda )=X_2.$$
(61)
Differentiating these, we see that

$$X_1B(\lambda )=X_2B(\lambda )=0,$$
(62)
and these also hold for derivatives of $`B`$ with respect to $`\lambda `$.

On differentiating (47) we find that (dropping the subscripts on $`R`$ and $`Y`$)

$$X^{}=R(YB+X_2nX_1X_1nX_2+(12\lambda )X_1X_2n)\stackrel{~}{R}.$$
(63)
It follows that

$$X^{}|_{\lambda =0}=X_2nX_1X_1nX_2+X_1X_2R_0n\stackrel{~}{R}_0,$$
(64)
where $`R_0=R(\lambda =0)`$. The tangent vector at $`X_1`$ is therefore given by (up to a scale factor)

$`T_1`$ $`=X_1nX_2X_1nX_1X_2(R_0n\stackrel{~}{R}_0)X_1n`$
$`=\left(R_0\left((X_1X_2n)X_1\right)\stackrel{~}{R}_0\right)n`$
$`=(C_1X_1)n.`$ (65)
But this is precisely the tangent vector to the circle $`C_1`$ at that point $`X_1`$. So the tangent at each of the control points is defined by the circle through the control point and its two adjoining points. The tangent vectors at a connection are therefore continuous, which ensures $`G_1`$ continuity.

Next we need to consider $`G_2`$ continuity. For this we need the circle

$`C_v`$ $`=XX^{}X^{\prime \prime }`$
$`=R\left(Y(Y^{}+YB)(Y^{\prime \prime }+2Y^{}B+YB^{})\right)\stackrel{~}{R}.`$ (66)
Evaluating this at $`\lambda =0`$ we obtain

$`C_v`$ $`=R_0(X_2nX_1)(X_2nX_1X_1nX_2+X_1X_2n)`$
$`(2X_1X_2n+2X_1X_2nB(0))\stackrel{~}{R}_0`$
$`=2(X_1nX_2nX_1X_2)C_1`$
$`2(X_2nX_1X_2)R_0X_1(X_1nX_2+X_1X_2n)(nB(0))\stackrel{~}{R}_0.`$ (67)
This contains a term in $`C_1`$, and an additional term controlled by the bivector $`B`$ (evaluated at $`\lambda =0`$). The first term is the desired one, as it ensures that the radius of curvature at the control point is defined by the circle through it. The second term is not wanted, and implies that the simple scheme defined by equation (42) does not guarantee $`G_2`$ continuity. This was first pointed out by Séquin & Lee . To provide $`G_2`$ continuity we need only ensure that the derivative of the interpolating rotor $`R`$ vanishes at the endpoints. This is simply achieved by replacing equation (42) with

$$\widehat{C}_{12}(\lambda )=\frac{1}{\mathrm{sin}(\theta )}\left(\mathrm{sin}\left((13\lambda ^2+2\lambda ^3)\theta \right)\widehat{C}_1+\mathrm{sin}((3\lambda ^22\lambda ^3)\theta )\widehat{C}_2\right).$$
(68)
This blending scheme does now have $`G_2`$ continuity, and is the scheme use to produce the figures in section 5. In figure 8 we show a comparison between the two schemes. As expected, the $`G_2`$ scheme smoothes out some of the curvature around the control points. One can extend this idea to obtain any desired order of continuity. This is achieved by replacing the polynomial in equation (68) by a higher-order blend.

## 8 Sphere blending

One can straightforwardly apply the ideas developed here to swept surfaces using a version of the scheme described by Szilvási-Nagy & Vendel . But in this section we aim to explore an alternative idea, based on sphere blending. A sphere is described as a 4-vector in the conformal geometric algebra, and these can be transformed and interpolated in a similar manner to circles.

As a simple example, suppose that the points $`X_1`$, $`X_2`$ and $`X_3`$ define the vertices of a triangle. To each corner we attach a sphere, $`S_1`$, $`S_2`$ and $`S_3`$. Each sphere passes through the three vertices of the triangle, and a fourth point which can be viewed as a control point. That is, we can write

$`S_1`$ $`=A_1X_1X_2X_3`$
$`S_2`$ $`=A_2X_1X_2X_3`$
$`S_3`$ $`=A_3X_1X_2X_3.`$ (69)
These assume that none of $`A_1`$, $`A_2`$ and $`A_3`$ lie on the circle defined by $`X_1`$, $`X_2`$ and $`X_3`$. We next need to define a blend over the surface. For this we introduce the barycentric coordinates $`(\lambda ,\mu ,\nu )`$, subject to

$$0\lambda ,\mu ,\nu 1,\lambda +\mu +\nu =1.$$
(70)
We let $`Y(\lambda ,\mu )`$ denote the conformal representation of the point in the triangle corresponding to the barycentric coordinates $`(\lambda ,\mu ,1\lambda \mu )`$. Similarly, we define a (linear) sphere blend by

$$S(\lambda ,\mu )=\lambda \widehat{S}_1+\mu \widehat{S}_2+\nu \widehat{S}_3.$$
(71)
We could employ a trigonometric blending scheme over the (abstract) sphere defined by the unit 4-vectors $`\widehat{S}_1`$, $`\widehat{S}_2`$ and $`\widehat{S}_1`$, but this raises a number of complications. There is no single straightforward generalisation of barycentric coordinates over a spherical triangle, and each alternative scheme has its own merits and drawbacks . Here we have adopted the simplest, linear blending scheme, which generates interesting surfaces.

Now that we have the sphere and the point on the triangle defined, all that remains is to define the conformal transformation from the point to the blended surface. First we define the plane $`P`$ through the three base points,

$$P=X_1X_2X_3n.$$
(72)
Next we define the rotor $`R`$ for the conformal transformation between the plane and the blended sphere,

$$R(\lambda ,\mu )=1\widehat{S}(\lambda ,\mu )\widehat{P}.$$
(73)
Finally, the surface itself is defined by the points $`X(\lambda ,\mu )`$, where

$$X(\lambda ,\mu )=R(\lambda ,\mu )Y(\lambda ,\mu )\stackrel{~}{R}(\lambda ,\mu ).$$
(74)
A typical surface defined in this manner is shown in figure 9. The result is aesthetically quite pleasing, yielding a smooth blend free of sudden changes in curvature.

Much work remains to extend this idea to a complete framework for defining blends over surfaces. The control points need to be chosen so as to reflect the geometry around the triangle, and continuity between blended spheres may be hard to achieve. Here we hope to have demonstrated that such an approach may be feasible.

## 9 Summary

Smooth splines between control points can be defined in terms of blended circles. The splines pass through all points, and have no extra control points. The natural geometric framework for handling circles is conformal geometry. This geometry is encapsulated in a simple, unified manner by employing the geometric algebra of a five-dimensional space. All transformations are easily defined, and algorithms can be written in such a way as to minimise problems with special cases. Similar ideas can be applied to sphere blends over a surface, and in future work we will explore the potential of this idea for design and for encoding surface data.

## Acknowledgements

CD thanks the EPSRC for their support.


format: text

text:
# The Dwarf Spheroidal Companions to M31: Variable Stars in Andromeda II1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., (AURA), under NASA Contract NAS 5-26555.

## 1 Introduction

There have been many studies of the stellar populations of the Galactic dwarf spheroidal (dSph) galaxies (see Mateo 1998 and references therein) which have enabled astronomers to better understand chemical evolution and star formation histories in less complex environments than for spiral or elliptical galaxies. In particular, the nature of dwarf galaxies has become important for understanding the formation of galaxies in general, as it is now believed that the halos of more massive galaxies were formed, at least in part, by the “cannibalism” of dwarf galaxies. This process continues today, with the accretion of the Sagittarius dSph galaxy into the halo of our own Galaxy being a prime example (Ibata, Gilmore, & Irwin 1994).

An important method for investigating the properties of dSph galaxies is to determine the parameters of their variable stars. For example, the simple presence of RR Lyrae stars (RRLs) is indicative of an older stellar population (age $`>10`$ Gyr). Due to their nearly uniform luminosities, RRLs can also be used to determine the distance to the system in which they are found. Every dSph galaxy surveyed for variable stars has also been found to contain at least one anomalous Cepheid (AC). This type of Cepheid variable derives its name from having a period-luminosity relationship that does not follow either the classical Cepheid or Population II Cepheid relationships. ACs are believed to be either young, metal-poor stars or stars that have formed from mass-transfer in a binary system (see, for example, the discussion in Pritzl et al. 2002, hereafter Paper I).

The dSph companions of the Milky Way Galaxy have all been surveyed for variable stars. Recently, we have begun a variability survey of the dSph companions to the Andromeda Galaxy using the Hubble Space Telescope (HST). We found the dSph galaxy Andromeda VI (And VI) to contain a significant population of RRLs and a number of ACs (Paper I). In this paper we present the results of our survey of Andromeda II (And II). And II, while having a mean metallicity slightly higher than that of And VI, has been shown to have a substantial metallicity spread (Côté, Oke, & Cohen 1999; Da Costa et al. 2000, hereafter DACS00). Here we show that the RRLs in And II also appear to have a broad metallicity distribution. We also discuss some of the other properties traced out by the RRLs and present the discovery of one AC. The relation of the AC specific frequency and the mean metallicity of the parent galaxy is also examined.

## 2 Observations and Reductions

Observations of And II were taken by the HST Wide Field Planetary Camera 2 (WFPC2) as part of the GO Program 6514 on 1997 August 29 and 1997 September 3, using the same orientation. Three 1200 s integrations through the F555W filter and seven 1300 s integrations through the F450W filter were taken during the first set of observations. Identical exposure times were used for the second set of observations, but there were four F555W and eight F450W integrations. The images were offset from the center of the dSph galaxy in order to avoid bright foreground stars (cf. Fig. 1 of DACS00).

The point-spread function fitting photometry was performed with allframe (Stetson 1994) as described in Paper I. The aperture and charge-transfer efficiency corrections along with the $`B`$,$`V`$ calibrations were carried out in the same fashion as outlined in Paper I. The main difference in the calibration of the And II dataset is the availability of the photometry from the combined WFPC2 images by DACS00, which form the basis of their color-magnitude diagram (CMD). We were able to compare our photometric calibrations with that of DACS00.

The results of this comparison are presented in Table 1. Here we list, for each wide-field CCD, the mean magnitude difference (in the sense this study – DACS00) and the standard deviation of the mean of the differences, along with the magnitude range over which the comparison was made, and the number of stars. Candidate variables have been excluded from the calculation. There are obvious systematic offsets between the two sets of photometry but it is reassuring to note that there are no signs of any trend with magnitude. The origin of these differences is unclear although, given the different photometry techniques and calibration processes employed (see Paper I and DACS00), the existence of such zero point differences is not particularly surprising. DACS00 demonstrated that their WFPC2 photometry agrees well with independently calibrated ground-based photometry for stars in the field of And II. Consequently, we have chosen to adjust the photometry of this study by the mean differences given in Table 1 to align it with that of DACS00.

## 3 Variable Stars

As for Paper I, we searched the photometry data for variable stars using a routine created by Dr. Peter Stetson called daomaster. This routine compares the rms scatter in the magnitudes to that expected from the photometric errors returned by the allframe program. There were only a handful of RRLs found on the PC and no ACs. Because the PC covers only a small area of sky compared to the WF CCDs, and because of the lack of sufficiently bright stars for reliable aperture corrections on the PC, we do not include the PC data in our results.

Figure 1 shows the location of the variable stars in the CMD created by DACS00. We removed the DACS00 photometric data from the CMD for those stars which were found to be variable in this paper and instead plotted them as plusses representing the mean magnitude and color of the variable star. The RRLs can be seen filling in the gap in the CMD along the horizontal branch (HB). About 1.2 mag brighter than the HB is the single AC we found in this survey. The periods, magnitudes, and light curves were determined using the programs created by A. C. Layden as outlined in Paper I. We note that due to the method we have used to calibrate the data to the $`B`$, $`V`$ system, the colors and magnitudes for the variable stars have inherent uncertainties on the order of 0.02-0.03 mag. Not withstanding this fact, we find that our method allows us to analyze the general properties of the And II variable star population.

The stars that lie to the blue of the red giant branch and those DACS00 stars found within and near the instability strip on the HB in Figure 1 were examined for variability, but no additional variables were found. Two plausible explanations for the non-variability, especially for those stars found within the instability strip, are either (a) these stars are not variables, but errors in the photometry, particularly in the colors, have scattered them into an apparent location in the instability strip, or (b) these are indeed variables but the limited time coverage of the observations meant that we did not detect the variability.

DACS00 performed an initial variable star search of And II. In their paper they report detecting approximately 30 candidate variable stars. Fig. 3 of DACS00 shows the location of these variables in the And II CMD and it shows a lot of dispersion among the RRLs. As discussed in Paper I the reason for such a dispersion is that the magnitudes and colors plotted by DACS00 come from 2 pairs of combined frames which, consequently, do not correctly sample the variable star light curves. They also show light curves for four of the RRLs they detected (see Fig. 6 of DACS00). These RRLs match V57, V60, V05, and V33 in this survey, going from the top-down in their figure. The periods for V57 and V05 match well in both studies, while the periods found for V60 and V33 in this survey and DACS00 are aliases of each other. The method used in this paper to obtain the photometry for the variable stars results in having more data points than were available to DACS00 allowing a more reliable choice of the period here.

A total of 73 variable stars were found on the WFC chips, one of which is an AC. This star will be discussed in the following section. Table 2 presents the photometric properties for these variable stars while their individual photometric $`B`$ and $`V`$ data are in Tables 3 and 4. Column 1 of Table 2 lists the star’s ID, while the next two columns give the RA and Dec. Finding charts for the variable stars are shown in Figure 2. Column 4 lists the period of each star. The intensity-weighted $`V`$ and $`B`$ magnitudes along with the magnitude-weighted colors $`(\mathrm{B V})_{\mathrm{mag}}`$ are shown in columns 5-7. Columns 8 and 9 give the $`V`$ and $`B`$ amplitudes of the variable stars. The remaining columns will be discussed later in the paper. In Figure 3 we present the light curves for all of the variable stars.

## 4 And II Anomalous Cepheid

We were able to find one AC, V14, in the field-of-view of the WFPC2. This follows the trend of at least one AC being found in every dSph galaxy surveyed for variable stars. ACs are believed to be either stars that have increased mass due to mass transfer in a binary system (Renzini, Mengel, & Sweigart 1977) or stars from an intermediate age population (Demarque & Hirshfeld 1975; Norris & Zinn 1975) although metal-poor abundances are required in both cases. The main way to recognize ACs is by their higher luminosity when compared to the RRLs in a system, typically 0.5-2.0 mag brighter than the HB.

The ACs follow a period-luminosity relation different than that for classical Cepheids or Population II Cepheids, hence the name “anomalous.” In Paper I we revised the period-luminosity relations for the ACs (see §4.1 in Paper I). Using the data from Table 2 and given $`(mM)_0=24.17\pm 0.06`$ and $`E(\mathrm{B V})=0.06\pm 0.01`$ for And II from DACS00, we find $`A_V=3.1E(\mathrm{B V})=0.19\pm 0.03`$ resulting in $`M_V=0.75`$ mag and $`M_B=0.43`$ mag for V14. This places the AC in And II with the other known ACs (see Table 4 of Paper I) in the period-luminosity diagram in Figure 4. V14 falls nicely along the first-overtone AC line in both $`M_B`$ and $`M_V`$.

We note that the light curve shape of V14 is quite asymmetric which would lead one to think it is pulsating in the fundamental mode. Yet, it is clearly a first-overtone mode pulsator as seen in Figure 4. This reinforces the view we expressed in Paper I that the shape of the light curve isn’t as clear an indicator of pulsation mode for an AC as is its location in a plot of the absolute magnitude versus the logarithm of the period. Given the $`B`$ and $`V`$ amplitudes for V14 in Table 2, this star also falls nicely among the other first-overtone mode AC pulsators in the period-amplitude diagram shown as Fig. 7 of Paper I.

## 5 And II RR Lyrae

RRLs are useful when examining the stellar populations of a system. Their periods, amplitudes, and magnitudes are probes of the distance, metallicity, and age of the system in which they belong. We found 72 RRLs in our field-of-view of And II, with 64 pulsating in the fundamental mode (RRab) and eight pulsating in the first-overtone mode (RRc). Their mere presence indicates that there are stellar populations of age $`>10`$ Gyr in And II. Due to the low amplitudes of the RRc stars and the uncertainties in our photometry, it is likely that we did not detect all of this type of RRLs. We compare the period distribution of the RRL population of And II to those of other dSph galaxies in Figure 5 using the sources listed in Table 6 of Paper I and the values derived above for And II. The mean metallicity increases from the top, Ursa Minor, to the bottom, Fornax. It should be noted that although the mean metallicity of Fornax is $`[\mathrm{Fe}/\mathrm{H}]=1.3`$ dex (Mateo 1998), the RRLs may originate from the metal-poor end of Fornax’s metallicity distribution (Bersier & Wood 2002). Looking at the RRab stars, there appears to be an overall trend of these stars shifting toward shorter periods as the metallicity increases. And II shows more of a spread in period among its RRab stars as compared to the other dSph galaxies. It is uncertain how much of this distribution may be due to the spread in metallicity that is known to exist in this dSph galaxy ($`\sigma _{\mathrm{int},[\mathrm{Fe}/\mathrm{H}]}0.36`$ dex, DACS00).

The mean periods for the RRLs are $`P_{ab}=0.571`$ day and $`P_c=0.363`$ day. In comparison to the period-metallicity relationship defined by the Galactic globular clusters, the mean period for the RRab stars is consistent with the mean metallicity in And II as determined from the red giant branch. For the RRc stars, however, the mean period is slightly longer than that expected if these stars have the same mean metallicity as the red giants. As discussed above, the RRc mean period may be unreliable as it is possible we are missing many of these stars. Nevertheless, the ratio of RRc stars to the total number of RRLs, $`N_c/N_{RR}`$, is 0.11, consistent with the relatively high mean metallicity of And II. If we place And II in Table 6 of Paper I using the mean metallicity from the red giant branch, it follows the trend in which the mean period for the RRab stars decreases as the metallicity increases. To illustrate this, we reproduce the metallicity versus mean RRab period diagram of Paper I (see Fig. 10 of that paper) in Figure 6, but now include And II. With the exception of $`\omega `$ Centauri and M2, there is a clear trend of increasing mean period with decreasing $`[\mathrm{Fe}/\mathrm{H}]`$ for the Galactic globular clusters. As noted in Paper I, the mean periods of the RRab stars in the dSph galaxies fill in the gap between the metal-rich and metal-poor Galactic globular clusters. And II clearly also follows this trend – it lies in Figure 6 among the Galactic globular clusters which have similar metallicities to the And II mean abundance.

We find the mean magnitude of the And II RRLs to be $`V=24.87\pm 0.03`$ mag, where the uncertainty is the aperture correction uncertainty, the photometry zeropoint uncertainty, and the spline-fitting uncertainty added in quadrature to the standard error of the mean. This matches the estimate of the level of the HB by DACS00, $`V_{\mathrm{HB}}=24.93\pm 0.03`$ mag based on the mean magnitude of the HB stars. We used the same equation as DACS00 to estimate the absolute magnitude of the RRLs,

$$M_{V,\mathrm{RR}}=0.17[\mathrm{Fe}/\mathrm{H}]+0.82$$
(1)
from Lee, Demarque, & Zinn (1990). We find $`M_{V,\mathrm{RR}}=+0.57`$ and $`A_V=3.1E(\mathrm{B V})=0.19\pm 0.03`$ given $`[\mathrm{Fe}/\mathrm{H}]=1.49\pm 0.11`$ and $`E(\mathrm{B V})=0.06\pm 0.01`$ from DACS00. As a result, we estimate the distance to And II to be $`665\pm 20`$ kpc, which matches the estimate made by DACS00 of $`680\pm 20`$ kpc.

### 5.1 RR Lyrae Period-Amplitude Diagram

The period-amplitude diagram provides an excellent diagnostic of the properties of the RRLs since neither of these quantities is dependent on distance or reddening. It is generally true that the position of RRLs in a period-amplitude diagram is dependent on its metallicity (Sandage 1993b). Still, there are other factors besides metallicity, such as evolution and age, that may shift a RRL’s position in the diagram. There is a clear division between the shorter period, smaller amplitude RRc stars and the RRab stars in the And II period-amplitude diagram (Figure 7). A wide spread in $`\mathrm{log}P`$ is seen among the RRab stars for a given $`A_B`$, that may or may not be due to the spread in $`[\mathrm{Fe}/\mathrm{H}]`$ for And II (DACS00).

There may also be some aliasing in the periods. As with And VI (Paper I), we used the period-finding routines created by Dr. Andrew Layden (Layden & Sarajedini 2000 and references therein) to find the best period for each variable. We also asked Dr. Gisella Clementini to determine what periods she would find for a small number of our variable stars using her period-finding program GRATIS (GRaphical Analyzer of TIme Series, see Clementini et al. 2000, 2003). The periods she found matched those determined from Andrew Layden’s program to $`\pm 0.001`$ day. For And VI we were able to make use of the period-amplitude diagram to see if the period adopted for each star was reasonable. However, the red giant branch population of And II shows a large abundance dispersion which may well also occur among the RRLs. This limits the usefulness of the period-amplitude diagram in determining the most appropriate period. Therefore, we have adopted for each star the best-fit period that was returned from Layden’s period-fitting program with no further screening.

Sandage, Katem, & Sandage (1981) detected a shift in period between the RRLs in M3 and M15 in the period-amplitude diagram. Sandage (1981; 1982a,b) later found that this period shift was dependent on the metallicity of a cluster where metal-poor clusters tend to have longer periods for a given amplitude when compared to metal-rich clusters. Using M3 as the fiducial cluster with zero period shift, Sandage (1982a,b) found

$$\mathrm{\Delta }\mathrm{log}P=(0.129A_B+0.088+\mathrm{log}P).$$
(2)
Sandage found that more metal-rich clusters have $`\mathrm{\Delta }\mathrm{log}P0.01`$, while the more metal-poor clusters have $`\mathrm{\Delta }\mathrm{log}P0.05`$. Between these two values, a gap exists where few RRLs are found. On the other hand, as discussed in Paper I (see §5.2), dSph galaxies have been shown to contain RRLs within the $`0.05\mathrm{\Delta }\mathrm{log}P0.01`$ range (Sextans: Mateo et al. 1995; Leo II: Siegel & Majewski 2000; Andromeda VI: Paper I). We list in column 10 of Table 2 the $`\mathrm{\Delta }\mathrm{log}P`$ values for the RRab stars of And II. In Figure 8 we plot these values against $`\mathrm{log}P`$ for each RRab star. It is clear in the figure that a number of And II RRab stars are found in-between $`0.05\mathrm{\Delta }\mathrm{log}P0.01`$. We find that a majority of the RRab stars are found in the metal-rich region of the diagram ($`\mathrm{\Delta }(\mathrm{log}P)>0.01`$). This is consistent with the mean metallicity of And II. When comparing And II to And VI (see Fig. 12 of Paper I), we see that there is a larger range in $`\mathrm{\Delta }\mathrm{log}P`$ values in And II than in And VI. This suggests that there is a sizable abundance spread among the And II RRLs, a result that is perhaps not unexpected given the large abundance dispersion among the red giants in this dSph galaxy. The mean period shift for And II is 0.01, which means that the RRab in And II are slightly more metal-rich than M3. Using

$$\mathrm{\Delta }\mathrm{log}P=0.116[\mathrm{Fe}/\mathrm{H}]+0.173$$
(3)
from Sandage (1982a), we find $`[\mathrm{Fe}/\mathrm{H}]=1.41`$ for And II. This estimate of the mean metal abundance of the And II RRLs is consistent with the mean abundance for And II red giants found by DACS00, $`[\mathrm{Fe}/\mathrm{H}]=1.49`$.

### 5.2 Metallicity Estimates from RR Lyrae

Alcock et al. (2000) developed a relation between the metallicity of a RRab star and its period and $`V`$-band amplitude. This was done in order to estimate the metallicity of the RRab stars in the Large Magellanic Cloud. Using the RRab stars in M3, M5, and M15, Alcock et al. found,

$$[\mathrm{Fe}/\mathrm{H}]_{\mathrm{ZW}}=8.85(\mathrm{log}P_{ab}+0.15A_V)2.60,$$
(4)
where ZW refers to the Zinn & West (1984) scale. The metallicity of the RRab stars in those systems was predicted with an accuracy of $`\sigma _{[\mathrm{Fe}/\mathrm{H}]}=0.31`$ per star. For the LMC, Alcock et al. found good agreement between the median metallicity from their relation and those from previous estimates. We applied this relation to our And VI RRLs (see §5.3 of Paper I) and found the mean of the RRL metallicity distribution generated from this relation to agree well with the mean metallicity for the red giants determined by Armandroff, Jacoby, & Davies (1999) using the red giant branch mean $`VI`$ color. The width of the metallicity distribution, however, was essentially identical to that expected from the uncertainty ($`0.3`$ dex) in the individual abundance determinations using this relation. Consequently, from these data, there is no evidence to support (or rule out) an abundance spread among the And VI RRLs.

We applied this equation to the RRab in And II and list the individual metallicities in column 11 in Table 2. Figure 9a shows the distribution of the individual $`[\mathrm{Fe}/\mathrm{H}]`$ from the RRab in And II. There is a wide range among the metallicity values similar to the histogram of metallicity values derived from red giant colors for And II as found by DACS00 and shown in Figure 9b. We also plot in Figure 9c the \[Fe/H\] distribution from the RRab stars in And VI (Paper 1). This figure emphasizes the differences between the metallicity distributions. Where the width of the \[Fe/H\] distribution for the And VI RRab stars was on the order of the uncertainty of the Alcock et al. equation ($`\sigma 0.3`$), the width for the And II stars is clearly larger than this.

Due to the wide spread in the metallicities, we did not attempt to fit a gaussian to the distribution as was done for And VI in Paper I. However, we find the mean to be $`[\mathrm{Fe}/\mathrm{H}]=1.46`$ and the median to be $`[\mathrm{Fe}/\mathrm{H}]=1.36`$. These both agree well with the mean metallicity found by DACS00 using the colors of the red giant branch stars, $`[\mathrm{Fe}/\mathrm{H}]=1.49\pm 0.11`$. It appears that the mean metallicity values resulting from RRLs and red giants in And II are consistent. Even though the metallicity distributions for the RGB and RRab stars are comparable, there is no reason why one should expect this. The mapping from the RGB to the instability strip is a complex function of (at least) metallicity, age distributions, and evolution within the instability strip as well as whatever parameters control mass-loss. Most And II RGB stars do not become a RRL, which is clear from looking at the stellar distribution along the HB. Nevertheless, it is intriguing that the range of abundances from the RGB and RRab stars is similar and they both show a peak in the distribution near the metal-rich end.

Studying cluster and field RRLs with a wide range of metallicities, Sandage (1993a) related the metallicities to the average periods of RRLs through the relations

$$[\mathrm{Fe}/\mathrm{H}]_{\mathrm{ZW}}=(\mathrm{log}P_{ab}0.389)/0.092$$
(5)
$$[\mathrm{Fe}/\mathrm{H}]_{\mathrm{ZW}}=(\mathrm{log}P_c0.670)/0.119.$$
(6)
Siegel & Majewski (2000) and Cseresnjes (2001) found that these relations applied well to the RRLs in dSph galaxies. We found in Paper I that the RRab estimate for And VI was within the combined errors of the estimate and that of Armandroff, Jacoby, & Davies (1999). The RRc abundance estimate, however, was lower than expected from the RRab and red giant branch stars, perhaps as a consequence of incompleteness in the sample of RRc stars, especially at shorter periods. The mean periods for the RRLs in And II are $`P_{ab}=0.571\pm 0.005`$ day and $`P_c=0.363\pm 0.005`$ day giving $`[\mathrm{Fe}/\mathrm{H}]_{\mathrm{ZW}}=1.58\pm 0.04`$ (internal error) for the RRab and $`[\mathrm{Fe}/\mathrm{H}]_{\mathrm{ZW}}=1.93\pm 0.05`$ (internal error) for the RRc. The metallicity estimate from the RRab stars, although slightly lower than the estimate made by DACS00, is within the combined errors. Once again, however, the estimate from the RRc stars is lower than that for the RRab and the red giant branch stars. It is possible that we have not detected all the RRc stars, especially at shorter periods, resulting in an inappropriate mean period (see §5).

## 6 The Frequency of Anomalous Cepheids in dSph Galaxies

Mateo, Fischer, & Krzeminski (1995) examined the specific frequency of ACs (the number of ACs per $`10^5L_{V,}`$) in the Galactic dSph galaxies and found a strong correlation between the specific frequency of ACs and the total luminosities of the parent galaxies: low luminosity dSph galaxies have higher specific frequencies of ACs. Interpretation of this correlation, however, is complicated. For example, Draco, Ursa Minor, and Carina all have similar luminosities and similar AC specific frequencies, but they have very different stellar populations. Draco and Ursa Minor are dominated by old stars and thus, in these particular dSphs, the ACs must have their origin in the binary mass transfer mechanism. On the other hand, Carina has a strong intermediate-age population and consequently, at least some of the ACs in this dSph may originate as single younger stars, rather than from (old) mass transfer binaries. Given these differences, there is no obvious reason why Carina should end up with a similar specific frequency to Draco and Ursa Minor. Similarly, Leo I has a dominant and relatively young intermediate-age metal-poor population (Gallart et al. 1999), yet its specific frequency of ACs is notably lower than that of Draco, Ursa Minor and Carina. Mateo et al. (1995) noted that dSph luminosities are strongly correlated with mean abundance (lower abundances for lower luminosities) and that this might be the underlying cause of the AC specific frequency - luminosity correlation: regardless of the origin of the ACs, they occur more frequently in metal-poorer environments.

There are now several more surveys for variable stars in dSph galaxies, including our surveys of the M31 dSph galaxies, than were available to Mateo et al. (1995). This may help disentangle the independent parameter(s) that drive the AC frequencies. In Table 5, we present revised values for Mateo et al. Table 7; we include the results from the latest surveys. For example, Dall’Ora et al. (2003) have identified 15 ACs in Carina, almost twice as many as previously known. In Table 5, Column 1 lists the dSph galaxy and column 2 is the mean \[Fe/H\]. The absolute magnitude of the galaxy is in column 3, and the structural parameters for each galaxy are given in columns 4–6 as the position angle measured from the North to the East, the ratio of the semi-minor axis to the semi-major axis, and the exponential scale length. The latter parameters are taken from the references given in the notes to the table. The percentage of the dSph luminosity surveyed and the number of ACs found are listed in columns 7 and 8, respectively. The specific frequencies of the ACs are then given in column 9, along with the logarithms of the specific frequencies and their uncertainties in column 10. $`S`$ is defined as the number of ACs per $`10^5L_{V,}`$ and was calculated for only the fraction of the dSph’s luminosity that was surveyed for variable stars. The uncertainties in the $`\mathrm{log}S`$ values were determined by assuming that the number of observed ACs is subject to Poisson statistics.

In order to calculate the luminosity surveyed in And VI and And II, we first determined the orientation and location of the region covered by the WFPC2 images relative to the galaxy. We then calculated the intensity of the WFPC2 overlap region relative to the entire galaxy by approximating the surface brightness distribution of the galaxy by the relation $`I=I_0e^{(r/r_0)^n}`$ (Sersic 1968), and making the appropriate allowance for the ellipticities. The $`V`$-band surface brightness profile parameters were taken from Caldwell et al. (1992) for And II and from Caldwell (1999) for And VI. From the luminosity fractions just calculated and the total visual magnitudes of the galaxies (column 3 of Table 5) we then derive the ($`V`$-band) luminosity encompassed by our images of each galaxy to which the AC surveys apply. The same process was applied for the Fornax dSph galaxy according to the survey area of the Bersier & Wood (2002) survey.

Using the data from Table 5, the logarithm of the specific frequency of ACs is plotted against $`M_V`$ for each galaxy in Fig. 10a, while in Fig. 10b, the logarithm of the specific frequency is plotted against the mean \[Fe/H\] for each dwarf. The Galactic dSph galaxies are shown as open circles and the M31 dSph galaxies as filled circles. In each plot the solid line is an unweighted least squares fit to the Galactic dSph data only. We confirm in Fig. 10a that a $`\mathrm{log}S`$ – (visual) luminosity correlation exists (correlation coefficient is 0.86) for the Galactic dSph galaxies. The correlation between $`\mathrm{log}S`$ and mean abundance shown in Fig. 10b is not as strong (correlation coefficient 0.74) but it is nevertheless significant – a correlation test shows that there is only a 4% probability (2.3-sigma) that the distribution of $`\mathrm{log}S`$ vs \[Fe/H\] arises by chance. This supports the contention of Mateo et al. (1995) that the mean abundance of dSph galaxies may be the underlying factor that governs the frequency of occurrence of AC variables, at least for the Galactic dSphs.

The addition of data for the M31 dSphs then allows us to investigate, at least in a preliminary way, whether these correlations apply also to the AC populations in dSph galaxies other than those of the Galaxy. The points for And II and And VI in Figs. 10a and 10b, indicate that these galaxies do appear to follow the relations as defined by the Galactic dSphs, but results for additional M31 dSphs are required to strengthen this conclusion. We will return to this question in our third paper in this series, where we will explore the variable star population, including ACs, of the M31 companions And I and And III.

## 7 Summary & Conclusions

Although it is one of the more metal-rich dSph galaxies, the properties of the variable stars in And II agree well with the trends defined by other dSph galaxies. We have found 73 variable stars in And II using HST/WFPC2 observations, one of which is an AC. The AC is found to have a period, absolute magnitude, and amplitude consistent with other ACs. The 72 RRLs were made up of 64 RRab stars and 8 RRc stars with mean periods of $`P_{ab}=0.571\pm 0.005`$ day and $`P_c=0.382\pm 0.005`$ day. The mean $`V`$ magnitude of the RRLs was found to be $`24.87\pm 0.03`$ mag resulting in a distance of $`665\pm 20`$ kpc for And II on the Lee, Demarque, & Zinn (1990) distance scale.

Using relations from Sandage (1982a, 1993a) and Alcock et al. (2000), we find that the properties of the RRab stars in And II yield estimates of the mean abundance of these stars that are in good agreement with the mean metallicity determined from both the colors (DACS00) and the line strengths (Côté et al. 1999) of And II’s red giants. With these results And II then follows the general trend seen in the Galactic globular clusters and Galactic dSphs in which the mean RRab period decreases with increasing metallicity. In addition, the And II RRLs show considerable spread in the period-amplitude diagram. We take this as an indication that there is a sizable abundance spread among these stars and have attempted to quantify this spread through use of the Alcock et al. (2000) relation. This yields an abundance distribution for the RRL which is considerably broader than that expected from the errors in the individual abundance determinations and contrasts with the case for And VI (Paper I) where there was no evidence for any intrinsic abundance spread among the RRL stars. The And II RRL abundance distribution is approximately uniform and has some features in common with the large metallicity spread derived by Côté et al. (1999) and DACS00 for the red giant branch stars (see Fig. 9).

We have investigated the existence of a relationship between the specific frequency of ACs and mean metallicities for dSph galaxies. As originally suggested by Mateo et al. (1995), we find for the Galactic dSphs there is a clear trend for higher specific frequencies at lower abundances. The M31 dSph galaxies And II and And VI also appear to follow this trend. However, more information on the frequency of occurrence of ACs in additional M31 dSph galaxies is required before we can fully compare the two sets of dSphs in this parameter plane.

This research was supported in part by NASA through grant number GO-08272 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. We would like to thank Peter Stetson for graciously sharing his PSFs for the WFPC2 and for the use of his data reduction programs. Thanks to Andrew Layden for the use of his light curve analysis programs. We would also like to thank Gisella Clementini for testing our period determinations. Thank you to Mario Mateo for his insightful referee comments and suggestions.


text:
# The anomalous behavior of coefficient of normal restitution in the oblique impact

thanks: Present address: Department of Chemistry, Kyoto University, Sakyo-ku, Kyoto, Japan, 606-8502

## Abstract

The coefficient of normal restitution in an oblique impact is theoretically studied. Using a two-dimensional lattice model for an elastic disk and an elastic wall, we demonstrate that the coefficient of normal restitution can exceed unity and has a peak against the incident angle in our simulation. We also explain this behavior based upon a phenomenological theory.

The coefficient of normal restitution $`e`$ is introduced to determine the normal component of the post-collisional velocity in the collision of two materials. The coefficient $`e`$ is defined by

$$𝐯(\tau _c)𝐧=e𝐯(0)𝐧,$$
(1)
where $`𝐯(\tau )`$ is the relative velocity of the centers of mass of two colliding materials at time $`\tau `$ measured from the initial contact, $`\tau _c`$ is the duration of a collision, and $`𝐧`$ is the unit vector normal to the contact plane. Though some text books of elementary physics state that $`e`$ is a material constant, many experiments and simulations show that $`e`$ decreases with increasing impact velocityStronge (2000). The dependence of $`e`$ on the low impact velocity is theoretically treated by the quasi-static theory Kuwabara and Kono (1987); Morgado and Oppenheim (1997); Brilliantov et al. (1996). We also recognize that $`e`$ can be less than unity for the normal impact without the introduction of any explicit dissipation, because the macroscopic inelasticities originate in the transfer of the energy from the translational mode to the internal modes such as the vibrationsMorgado and Oppenheim (1997); Gerl and Zippelius (1999); Hayakawa and Kuninaka (2002).

While $`e`$ has been believed to be less than unity in most situations, it is recently reported that $`e`$ can exceed unity in oblique impactsSmith and Liu (1992); Calsamiglia et al. (1997); Louge and Adams (2002). In particular, Louge and AdamsLouge and Adams (2002) observed oblique impacts of a hard aluminum oxide sphere on a thick elastoplastic polycarbonate plate and found that $`e`$ grows monotonically with the magnitude of the tangent of the incident angle $`\gamma `$. In their experiment, Young’s modulus of the plate is $`100`$ times smaller than that of the aluminum oxide sphere. They also suggested that $`e`$ can exceed unity for the most oblique impacts.

In this letter, we demonstrate that our two-dimensional simulation of the oblique impact based on Hamilton’s equation has yielded an increasing $`e`$ with $`\mathrm{tan}\gamma `$ and $`e`$ exceeds unity at the critical incident angle. Finally, we explain our results by our phenomenological theory.

Let us introduce our numerical modelKuninaka and Hayakawa (2003). Our model consists of an elastic disk and an elastic wall (Fig. 1). The width and the height of the wall are $`8R`$ and $`2R`$, respectively, where $`R`$ is the mean radius of the undeformed disk. The both side ends and the bottom of the wall are fixed. We place $`800`$ mass points at random in a disk with the radius $`R`$ and $`4000`$ mass points at random in a wall for the disk and the wall, respectively. We connect each mass point with its neighbor mass points by the Delaunay triangulation algorithmSugihara (2001), and undeformed nonlinear springs are placed on all the connections.

Each mass point $`i`$ on the lower half boundary of the disk feels the force, $`𝐅(l_s^{(i)})=aV_0\mathrm{exp}(al_s^{(i)})𝐧_s^{(i)}`$, where $`l_s^{(i)}`$ is the distance between $`i`$-th surface mass point of the disk and the nearest surface spring of the wall, $`a=300/R`$, $`V_0=amc^2R/2`$, $`m`$ is the mass of each mass point, $`c`$ is the one-dimensional velocity of sound, and $`𝐧_s^{(i)}`$ is the unit vector normal to the connection between two surface mass points of the wall Kuninaka and Hayakawa (2003). We should note that the strong repulsion $`𝐅(l_s^{(i)})`$ is introduced to inhibit the penetration of the disk to the surface of the wallGerl and Zippelius (1999). Thus, the dynamical equation of motion for each mass point $`i`$ of the lower half boundary of the disk is described by

$$m\frac{d^2𝐫_i}{d\tau ^2}=\underset{j=1}{\overset{N_i}{}}\left\{k_a𝐱_{ij}k_b𝐱_{ij}^3\right\}+\mathrm{\Theta }(l_{th}l_s^{(i)})aV_0\mathrm{exp}(al_s^{(i)})𝐧_s^{(i)},$$
(2)
where $`𝐫_i`$ is the position of $`i`$-th mass point, $`N_i`$ is the number of mass points connected to $`i`$-th mass point, $`𝐱_{ij}`$ is the relative deformation of the spring between $`i`$-th and $`j`$-th connected mass points, $`k_a`$ and $`k_b=k_a\times 10^3/R^2`$ are the spring constants. Here we introduce the step function $`\mathrm{\Theta }(x)`$, i.e. $`\mathrm{\Theta }(x)=1`$ for $`x0`$ and $`\mathrm{\Theta }(x)=0`$ for $`x<0`$, and the threshold length $`l_{th}`$ which is the average of the natural lengths of the springs of the disk. For internal mass points, the last term of the right hand side of eq.(2) is omitted. In most of our simulations, we adopt $`k_a=k_a^{(d)}=1.0\times mc^2/R^2`$ for the disk and $`k_a=k_a^{(w)}=1.0\times 10^2mc^2/R^2`$ for the wall. We do not introduce any dissipative mechanism in this model. Thus, during a collision, a part of initial translational energy of the disk is distributed into the vibrational energy of the disk and the wall. It should be noted that the macroscopic dissipation can be interpreted as the irreversible transfer of the energy from the translational motion to the internal vibration. When we introduce explicit dissipations in the model and add the gravity to the disk, we have confirmed that the compression of the disk can be described by two-dimensional Hertzian contact theory Kuninaka (2004).

In this model, the roughness of the surfaces is important to make the disk rotate after a collisionKuninaka and Hayakawa (2003). We modify the positions of the surface mass points of the flat wall and the smooth disk by using normal random numbers whose average and standard deviation are $`0`$ and $`3\times 10^2R`$, respectively.

Poisson’s ratio $`\nu `$ and Young’s modulus $`E`$ of this model can be evaluated by adding the external force to stretch the rectangle of random lattice numerically. We obtain Poisson’s ratio $`\nu =(7.50\pm 0.11)\times 10^2`$ and Young’s modulus $`E=(9.54\pm 0.231)\times 10^3mc^2/R^2`$, respectivelyKuninaka and Hayakawa (2003).

We solve the dynamical equation of motion (2) for each mass point with the initial speed $`|𝐯(0)|=0.1c`$ and the incident angle $`\gamma `$, and determine $`e`$ for each $`\gamma `$ according to eq.(1). All the results in this letter are obtained by averaging the results of 100 disks with different configurations of mass points. We use the fourth order symplectic integrator with the time step $`\mathrm{\Delta }\tau =10^3R/c`$.

Figure 2 is the relation between $`e`$ and $`\mathrm{\Psi }_1(𝐯(0)𝐭)/(𝐯(0)𝐧)=\mathrm{tan}\gamma `$, where $`𝐭`$ is the unit vector vertical to $`𝐧`$. The cross point is the mean value and the error bar is the standard deviation of 100 samples for each $`\gamma `$. This result shows that $`e`$ increases with increasing $`\mathrm{\Psi }_1`$ to exceed unity, and has a peak around $`\mathrm{\Psi }_1=6.0`$. The behavior of $`e`$ having the peak is contrast to that in the experiment by Louge and AdamsLouge and Adams (2002).

Here, let us explain our results. Louge and AdamsLouge and Adams (2002) suggest that their results can be explained by the rotation of the normal unit vector $`𝐧`$ arising from the local deformation of the wall’s surface. Thus, we aim to determine the angle of rotation of the unit vector $`\alpha `$ at each $`\gamma `$ from the theory of elasticity.

Figure 3 is the schematic figure of a hard disk moving from the left to the right on a wall, where $`l|x_bx_a|`$. Assuming that $`l`$ is small compared to $`R`$, the perimeter of the contact is approximated by

$$f(x)=(xx_c)^2/2Ry_c,$$
(3)
where $`(x_c,y_c)`$ is the lowest position of the disk. To calculate $`\mathrm{tan}\alpha (f(x_b)f(x_a))/l`$, we need to know the ratio of $`|x_cx_a|`$ to $`l`$. From the theory of elasticityGalin (1953), this ratio can be estimated as

$$\frac{x_cx_a}{l}=1\theta \text{with}\theta =\frac{1}{\pi }\mathrm{arctan}\frac{12\nu }{\mu (22\nu )},$$
(4)
where $`\mu `$ is the coefficient of the friction. We evaluate $`\mu `$ by $`\mu |𝐉𝐭|/|𝐉𝐧|`$ with $`𝐉=M(𝐯(\tau )𝐯(0))`$, where $`M`$ is the mass of the disk. The cross points in Fig.4 represent $`\mu `$ calculated from our simulation, where $`\mu `$ has a peak around $`\mathrm{\Psi }_1=3.0`$. From eqs.(3) and (4), $`\mathrm{tan}\alpha `$ is given by

$$\mathrm{tan}\alpha =\frac{2\theta 1}{22\theta }\frac{|x_cx_a|}{R},$$
(5)
where we evaluate $`|x_cx_a|=0.55R`$ from the maximum value of the compression of the disk. From eqs.(4) and (5), we obtain the relation between $`\mathrm{\Psi }_1`$ and $`\mathrm{tan}\alpha `$.

Next, let us calculate $`e`$ from the relation between $`\mathrm{\Psi }_1`$ and $`\mathrm{tan}\alpha `$. We introduce the rotated unit vectors, $`𝐧_\alpha `$ and $`𝐭_\alpha `$, as $`𝐧_\alpha =\mathrm{cos}\alpha 𝐧\mathrm{sin}\alpha 𝐭`$ and $`𝐭_\alpha =\mathrm{sin}\alpha 𝐧+\mathrm{cos}\alpha 𝐭`$, respectively. By introducing $`e_\alpha (𝐯(\tau _c)𝐧_\alpha )/(𝐯(0)𝐧_\alpha )`$, we can express $`e`$ in terms of $`e_\alpha `$ as

$$e=\frac{e_\alpha +\mathrm{\Psi }_2^\alpha \mathrm{tan}\alpha }{1\mathrm{\Psi }_1^\alpha \mathrm{tan}\alpha },$$
(6)
where $`\mathrm{\Psi }_1^\alpha =(𝐯(0)𝐭_\alpha )/(𝐯(0)𝐧_\alpha )`$ and $`\mathrm{\Psi }_2^\alpha =(𝐯(\tau _c)𝐭_\alpha )/(𝐯(0)𝐧_\alpha )`$. $`\mathrm{\Psi }_1^\alpha `$ also can be rewritten as

$$\mathrm{\Psi }_1^\alpha =(\mathrm{\Psi }_1\mathrm{tan}\alpha )/(1+\mathrm{\Psi }_1\mathrm{tan}\alpha ).$$
(7)
On the other hand, in the oblique impact of slipping disks, $`\mathrm{\Psi }_2^\alpha `$ is given by

$$\mathrm{\Psi }_2^\alpha =\mathrm{\Psi }_1^\alpha 3(1+e_\alpha )\mu _\alpha $$
(8)
in the two-dimensional situationWalton and Braun (1986). In eq.(8), $`\mu _\alpha `$, defined by $`\mu _\alpha =|𝐉𝐭_\alpha |/|𝐉𝐧_\alpha |`$, is given by

$$\mu _\alpha =\frac{\mu +\mathrm{tan}\alpha }{1\mu \mathrm{tan}\alpha }.$$
(9)
To draw the solid line in Fig.2, at first, we calculate $`\mu `$ and $`\mathrm{tan}\alpha `$ for each $`\mathrm{\Psi }_1`$. By choosing a fitting parameter $`e_\alpha =0.95`$ and substituting eqs.(7)-(9) into eq.(6) we obtain $`e`$ as a function of $`\alpha `$, $`\mu `$, and $`\mathrm{\Psi }_1`$. All points are interpolated with the cubic spline interpolation method to draw the solid curve. Such the theoretical description of $`e`$ is qualitatively consistent with our numerical result, though the theoretical value is a little smaller than the observed value.

Now, let us consider how $`\mu `$ depends on $`\mathrm{\Psi }_1`$ based on a phenomenological argument. For simplicity, we replace the roughness on the surface of the wall by a periodic array of asperities. When the disk hits one asperity, a fraction of the energy is absorbed in the wall. We assume that the tangential velocity $`v_t^{(i1)}`$ changes to $`v_t^{(i)}=(1\eta )v_t^{(i1)}`$ when the disk hits the $`i`$-th asperity. Thus, if the disk contacts $`N`$ asperities during the impact, the tangential speed at the release point becomes $`v_t(\tau _c)=v_t^{(N)}=(1\eta )^Nv_t(0)`$. Here we estimate the number of contacted asperities during the collision as $`N=\rho l_{sl}`$, where $`\rho `$ is the number of the asperities in a unit length on the surface and $`l_{sl}`$ is the length of sliding which can be evaluated as $`l_{sl}=v_t(0)\tau _c`$ with $`\tau _c=\pi (R/c)\sqrt{\mathrm{ln}\left(4c/v_n(0)\right)}`$Hayakawa and Kuninaka (2002). Thus, the tangential impulse $`J_t𝐉𝐭`$ is approximated by

$$J_t=M\{(1\eta )^{\rho l_{sl}}1\}v_t(0)\eta M\rho l_{sl}v_t(0),$$
(10)
where we assume small $`\eta `$. Now, we should answer the question how to determine $`\eta `$. When the impact velocity is large enough, $`\eta `$ becomes smaller because the asperities are broken down when the disk hits them. Therefore we may assume the form $`\eta =\eta _0/(1+\beta v_t^2(0)/v_n^2(0))=\eta _0/(1+\beta \mathrm{\Psi }_1^2)`$ with the introduction of the dimensionless parameters $`\eta _0`$ and $`\beta `$. Taking into account $`J_n𝐉𝐧=M(1+e)v_n(0)`$ and the definition of $`\mu =J_t/J_n`$ we obtain

$$\mu =\frac{\pi \eta _0\rho R}{1+e}\frac{\mathrm{\Psi }_1}{1+\beta \mathrm{\Psi }_1^2}\sqrt{\frac{\mathrm{\Psi }_1^2}{1+\mathrm{\Psi }_1^2}\mathrm{ln}\left(40\sqrt{1+\mathrm{\Psi }_1^2}\right)}\frac{|𝐯(0)|}{c}.$$
(11)
Here we use $`\mathrm{cos}\gamma =\sqrt{1/(1+\mathrm{\Psi }_1^2)}`$, $`\mathrm{sin}\gamma =\sqrt{\mathrm{\Psi }_1^2/(1+\mathrm{\Psi }_1^2)}`$. We adopt the numerical results for the value of $`e`$ at each $`\mathrm{\Psi }_1`$ (Fig.2) and $`|𝐯(0)|/c=0.1`$. The solid curve in Fig.4 is eq.(11) with fitting parameters $`\beta =0.21`$ and $`\eta _0\rho R=0.18`$, which reproduces our numerical result. We do not claim that our simple argument explains the experimental result because of two fitting parameters. However, we would like to emphasize that our picture captures the essence of the behavior of $`\mu `$ in the oblique impact.

Let us discuss our result. First, we emphasize that the novel phenomena of $`e`$ exceeding unity are obtained from the local deformation of the soft wall for the oblique impacts of a hard disk. When we simulate the impact between a disk and a hard wall, for $`k_a^{(w)}=10\times k_a^{(d)}`$, $`e`$ fluctuates around a constant to exceed unity abruptly around $`\mathrm{\Psi }_1=4.5`$Kuninaka (2004). This tendency resembles the experimental results by Calsamiglia et al.Calsamiglia et al. (1997). Thus, for smooth increase of $`e`$ to exceed unity, the wall should be softer than the disk. In addition, it is important to fix the initial kinetic energy of the disk. So far, we have confirmed that $`e`$ does not exceed unity when $`\mathrm{\Psi }_1`$ is controlled by changing $`v_t`$ with fixed $`v_n`$Kuninaka and Hayakawa (2003).

The second, the initial velocity of the disk and the local deformation of the wall are much larger than those in the experimental ones in ref.Louge and Adams (2002). They cause the most significant difference between our result and their resultLouge and Adams (2002). Because of the high speed impact in our case, there is a peak of $`e`$ for small $`\gamma `$. In fact, our simulation with $`|𝐯(0)|=0.01c`$ shows the shift of the peak for larger $`\gamma `$. Therefore, we expect that our model reproduces the result of ref.Louge and Adams (2002) for the low impact speed. In addition, we have carried out simulations with a disk of $`400`$ mass points and a wall of $`2000`$ mass points to check the effect of the system size. Although there is a slight difference between the results, the data are also well reproducible by our phenomenological theory.

The third, the local deformation of the wall also affects the relation between $`\mu `$ and $`\mathrm{\Psi }_1`$. In early studies, it has been shown that $`\mu `$ depends on the impact velocity Gorham and Kharaz (2000); Louge and Adams (2002). In our simulation, $`\mu `$ has a peak around $`\mathrm{\Psi }_1=3.0`$. This behavior is interpreted as that the asperities are flattened for large $`v_t`$. Equation (11) indicates that $`\mu `$ can increase with increasing $`\gamma `$ as in ref.Louge and Adams (2002) if we choose a suitable set of $`\beta `$ and $`\eta _0\rho `$. The difference between the results of their experiment and our simulation may be explained by the choice of these parameters.

The fourth, we adopt the static theory of elasticity to explain our numerical results in this letter. However, it is important to solve the time-dependent equation of the deformation of the wall’s surface to analyze the dynamics of impact phenomena. The dynamical analysis will be our future task.

In final, we indicate that the friction coefficient is derived from our Hamiltonian model. This friction comes from the irreversible energy transfer of the macroscopic translational motion to the internal vibration. The irreversibility is indeed related to the second law of thermodynamics if the number of the internal degrees of freedom is infinite. However, the irreversibility in our system which includes only $`10^3`$ mass points in a disk is incomplete as we can see that the Hertzian contact theory cannot be recovered without introduction of the explicit dissipation. We believe that the complete treatment of inelastic collision of macroscopic materials without introduction of dissipations will be a fundamental subject of nonequilibrium statistical mechanics.

In summary, we have carried out the two-dimensional simulation of the oblique impact of an elastic disk on an elastic wall. We have found that $`e`$ can exceed unity in the oblique impact, which is attributed to the local deformation of the wall. We have estimated the magnitude of the local defomation $`\alpha `$ based on the static theory of elasticity and derived the relation between $`e`$ and $`\mathrm{\Psi }_1`$ by taking into account the rotation of the normal unit vector of the wall’s surface. The relation between $`\mu `$ and $`\mathrm{\Psi }_1`$ is also related to the local deformation and is explained by the simple argument.

###### Acknowledgements.

We would like to thank Y. Tanaka and S. Nagahiro for their valuable comments. We would also like to thank M. Y. Louge for his detailed introduction of their experiment to us and some valuable comments. Parts of numerical computation in this work were carried out at Yukawa Institute Computer Facility. This study is partially supported by the Grant-in-Aid of Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (Grant No. 15540393) and the Grant-in-Aid for the 21st century COE ”Center for Diversity and Universality in Physics” from MEXT, Japan.