A mathematical framework for the analysis and comparison of contact detection methods for ellipses and ellipsoids

The purpose of this research is to provide a framework for the analysis and comparison of contact detection algorithms for pairs of ellipses and ellipsoids. This work focuses primarily on the category of algorithms that are the most computationally efficient and can produce estimates of the separation and the penetration distance between ellipses and ellipsoids. Specifically, only analytic representations of the ellipses and ellipsoids are considered and contact detection for moving pairs of ellipsoids is not treated. The first contribution is a mathematical framework for the study of these algorithms, most notably with existence and uniqueness proofs for classes of contact detection algorithms, formal descriptions of the asymptotics of pairs of ellipses in close contact (or overlap), and a global analysis of constraints on the normals. The framework highlights the key role played by the different definitions of contact found in the literature, independent of the numerical strategies deployed to estimate the separation/penetration distance. Specifically, it is shown that all the studied algorithms can be expressed as minimization problems, with or without non-binding constraints on the normal(s) at the contact point(s), and that the constraints can be used to identify the global minima among the critical points in the minimization problem. Another contribution of this research, based on the mathematical framework introduced, is a better classification of the known algorithms. These algorithms are compared on established test problems, and their strengths and weaknesses are highlighted and explained in terms of their classification. Furthermore, this research provides comparisons in speed and stability between the most efficient algorithms in each category over a large sample size of test problems. Among the other contributions, this research describes inexpensive but effective initial estimates of the contact to be used in iterative algorithms. Finally, the usefulness of the new framework is illustrated with the introduction of a fast algorithm combining some new and old ideas.


Introduction
This research describes a unifying framework for the analysis and comparison of contact detection algorithms (CDAs) for pairs of ellipses and ellipsoids 1 . The published CDAs for pairs of ellipses tend to propose incremental improvements over past methods, to offer few comparisons with significantly different algorithms and to rarely distinguish their underlying mathematical problem from their numeri-solutions, which are generated by a random but representative sampling from the space of configurations. Finally, this research presents a novel CDA which outperforms previous algorithms in computational time by a factor of at least three, thereby illustrating how these new mathematical notions can lead to significant progress on the numerics of contact detection. In summary, this research sheds new light on an old problem and provides many new directions for future research.
CDAs are an essential component of many complex simulations such as those of granular flows in geo-mechanics, but also of hopper and chute flows in the chemical and food industries, of failure and fracture of brittle materials, and even of the forces of ice flows on offshore structures. In the previous applications, CDAs play a key role in the computation of inter-particle forces and the use of ellipsoids rather than spheres or cylinders and increase the applicability and reliability of the simulations. These numerical applications are often variants of the discrete element method (DEM) [7], or of molecular dynamics simulations (with or without Coulomb forces) [1,10,22]. On the other hand, some applications require the same accuracy in the estimation of the separation/penetration distance between ellipses, but do not necessarily need to estimate forces, such as in robotics, computer graphics (CAD/CAM) [14], data uncertainty analysis, computer vision, or virtual reality applications, to name a few. In those cases, ellipses often serve as effective bounding boxes of geometrically complex objects for which collision detection and avoidance is crucial.
Given our intended applications, this study will omit discussing contact detection between particles not defined analytically by an ellipse and will also ignore contact detection for rapidly moving pairs of ellipses. In particular, we are not interested by the contact detection problem when the ellipse is approximated by segments of circles [25,28,37,49], by grid or polar representation of particles [24], by the four arc approximation of ellipses [28,49], as a polyhedral surface [17], or in the most general form as a combination of NURBS [30,48]. The analysis is also limited to soft particles, that is to pairs of ellipses for which an accurate estimate of the separation/penetration distance is required, usually to compute forces between particles. This excludes a large set of methods used in static molecular assemblies and based on the Perram-Wertheim theory, and its many successful extensions. Furthermore, when the particles are moving (and rotating) rapidly with respect to their traveling distance, then the techniques of Wang et al. [26,51] are more appropriate than those studied here. Nonetheless, similarities to these last two techniques are discussed, respectively, in Sects. 3.5 and 3.9. In this respect, this research will focus on highly accurate and computationally effective methods of contact detection for pairs of soft ellipses in the quasi-static regime.
Although the study of granular media should consider particles of almost any shape, it is clear that ellipses are the simplest convex models of particles which can adequately describe both kinetic and rotational energies. In fact, studies have shown that packed assemblies of ellipses exhibit behaviors absent from packed assemblies of disks [8]. Unfortunately, current CDAs for pairs of ellipses can incur a computational overhead that is more than an order of magnitude higher than between pairs of spheres [54], with CDAs between pairs of particles with even more complex geometries requiring at least two orders of magnitude more than spheres. This implies that there are limits to the scale that DEM/MD simulations can attain and of the conclusions that can be reached about the influence of particle geometry on macroscopic soil mechanics. In this sense, faster and more accurate CDAs for pairs of ellipses would help to scale simulations closer to realistic scales. Although it is conjectured that some of the results in this research will extend to hyperquadrics, the contact detection problem for pairs of particles with more general (non-convex) geometries might not fit into a similar framework.
The main idea of the CDA studied here is to find two points associated with the two particles, and from these determine a contact zone. One can then use those two points to define a contact point between elliptic particles, a normal direction, and the forces, of equal size but opposite direction. The accuracy and stability of the CDA therefore directly influence the behavior of the short-and long-term dynamics of granular assemblies. We emphasize here that the notion of contact point and of normal vector is not uniquely defined for two overlapping ellipses, in contrast to disks. In fact, as two disjoint ellipses approach each other, their relative velocities and angular motion play a role in determining the evolution of the contact region. In conclusion, our view is that it is more important to propose a well-defined and relevant notion of contact point and separation/penetration distance than to formulate the most general definition which might ultimately be impractical to compute.
The framework developed in this research leads to a natural classification of published CDAs into categories based on their chosen definition of contact, which we now summarize. The framework of Sect. 3 begins by motivating and defining three separate notions of contact, all expressed as minimization problems. The naïve approach for soft ellipses is to identify the points, if any, belonging to the intersection of the two ellipses. This method, referred to as the intersection algorithm (IA) in Sect. 5.1, was analyzed by Rothenburg et al. in [42] but is well known to be unstable and therefore to be avoided. The second most natural definition of contact, in the case when the pair is non-overlapping, begins by identifying the pair of closest points on the ellipses, formalized as the minimum distance pair (MDP) in Lemma 3. When the pair is overlapping, the solution to the minimization problem associated with the MDP is in fact the same as for the IA. For overlapping ellipses, it is therefore more natural to extend the definition of MDP by introducing an additional constraint on the normals, that was previously non-binding for disjoint ellipses. This extended notion of MDP serves as the basis for the closest co-normal algorithm (CCN) of Sect. 5.6, first studied by Wellmann et al. [53]. Unfortunately, this problem is quite difficult to solve since it leads to a large coupled nonlinear system of equations. For that reason, many researchers have developed algorithms that split into two separate smaller problems: for each ellipse, search for the point on that ellipse closest to the other ellipse. Typically, closeness is measured with respect to the induced norm of the associated ellipse and leads to the minimum potential pair (MPP) of Lemma 4. This approach serves as the basis to several algorithms [12,34,46], which differ only in the solution process but not in the underlying problem definition. This research also examines the common normal algorithm (CNA) of Lin et al. [31]. As far as we know, the only attempt at a survey of these algorithms has been a dedicated chapter in a monograph [32]. Finally, it should be noted that many of the published algorithms have been renamed in this survey because the original names often led to significant misunderstandings about their contact definition (IA, MDP, MPP) and their numerical formulation.
The contact theory of Perram-Wertheim [41,47] can only be used to detect overlaps but cannot predict separation/penetration distance. This theory only estimates a proxy of contact separation/penetration, not two points on the ellipses, and therefore cannot be used to accurately estimate the force between two ellipses. Hence, we will only indicate some of its similarities to the MDP and the MPP in Sect. 3.5, namely that their contact potentials are expressed as minimization problems and require global constraints on the fields of normals associated with the particles. It was first constructed by Vieillard-Baron [47] in 1972 for pairs of ellipses before being generalized by Perram and Wertheim [41] into an effective tool for the computation of thermodynamic properties of dense static assemblies.
The engineering and computer science literature contains several variants of the contact detection problem that will not be studied here, but which nonetheless warrant mention. One of these is the broad phase search of contact detection where, for each particle, a preliminary list of neighbors is identified before a narrow phase search is completed between the particle and its neighbors, e.g., the type of accurate contact detection we are concerned with. Given the high cost of a CDA for a pair of particles during the narrow phase, an efficient broad phase search can be essential, such as the parallel algorithm over GPUs proposed by NVIDIA [36]. We also remark that there is a large and growing literature on contact detection for swarms of robots and drones [50], with or without centralized control. For the promise of driverless vehicles to become a reality, image recognition and contact detection algorithms will need to become more accurate and reliable.
The mathematical framework for contact detection introduced in this research was deemed necessary in order to unify the current state of the art on the topic, as past contributions were often based on the researchers' intuition and sometimes lacked precise definitions, formal proofs and exhaustive benchmarks. From a mathematical point of view, this meant it was necessary to distinguish the underlying contact detection problem, with its own existence and uniqueness results for its solution, from the numerical contact detection algorithm, that is required in order to estimate that solution. The theory presented in Sect. 3 will deal only with the underlying mathematics of the problem, while Sect. 5 will distinguish algorithms based on their underlying contact detection problem and will detail the many different approximations that have been used to solve that problem.
The rigorous framework for CDA introduced in this paper is an attempt to formalize key concepts in contact detection. The first part of the analysis is to motivate the introduction of three separate notions of contact in the form of minimization problems: the intersection set (IS), the minimum distance pair (MDP), and the minimum potential pair (MPP). Existence and uniqueness theorems are given for these notions, which are, rather surprisingly, absent from the literature. In the case of the MPP, it is found that non-binding constraints on the normals exist and these are formalized in the notion of the co-gradient locus of Theorem 6, previously introduced by [46] as the line of common slope. The co-gradient locus, being an intrinsically defined object associated with each pair of ellipses, is then used to clarify the notion of two ellipses in near-perfect contact. Theorem 7 states that if an intrinsically numerical constraint, as introduced in Definition 6, is satisfied, then the ellipses are guaranteed to be in a configuration consistent with the use case in DEM, namely with at most two points of intersection. As far as the mathematical theory is concerned, not all results have been demonstrated in 3D, but the remaining elements of theory are easily formulated and will be highlighted throughout Sect. 3. Although we have attempted to present a rigorous and self-contained presentation of the framework, many open questions remain and it is our hope that this framework will help to clarify and drive forward future research on this basic question.
The framework goes beyond the theory and contains numerical examples and procedures to compare the different CDA. The comparison between CDAs is particularly difficult because they may have different problem definitions, rely on different numerical algorithms to solve the same intermediate problems, have different implementations, use different tolerance definitions and values, etc. This explains why the literature rarely contains comparisons between more than two CDAs. The approach taken is to generate a large set of pairs of ellipses, appropriately sampled from the space of configurations, to solve the contact detection problem for each CDA, and extract statistical information from the results. Given that we cannot predict beforehand all the configurations of pairs of ellipses which might highlight the weaknesses of different CDAs, it is important to sample the entire space of configurations of pairs of ellipses, notably the ratios of areas, the aspect ratios, the separation/penetration distance, the relative orientations, and the locations of the contact points along the perimeter of the ellipses. The statistical data that are extracted are the mean and variance of the computational time, the accuracy (because the exact MPP is known by construction), and the number of iterations of the underlying nonlinear solvers. These numerical experiments are carried out in Sect. 6, but the detailed algorithm used to generate random pairs of ellipses is detailed in a second paper by the authors [29].
A final contribution of this research is a novel CDA that incorporates some new and old elements and illustrates how the framework presented here can lead one to construct new algorithms. The new CDA, which we have named the Steered Geometric Potential Algorithm (S-GPA) for reasons that will become clear later, is an algorithm to compute MPP using the mapping of Sect. 4.2 and Newton's method. It includes both a novel technique to generate an initial guess of the contact points and a novel constraint ensuring that the algorithm converges iteratively to the unique MPP. Both of these new elements are of independent interest. This new CDA is only briefly explained in Sect. 5.4, but a complete description of the algorithm can be found in [29].
This review paper is organized as follows. We describe in Sect. 2 the mathematical representations of ellipses and ellipsoids. Section 3 introduces a detailed mathematical formulation of the contact detection problem, thereby laying the foundation results for our description of the various techniques in Sect. 5. Section 4.3 briefly reviews efficient initial estimates of the contact point. Afterward, we briefly review the methods and analyze their respective advantages and disadvantages in two-and three-dimensional problems. The last section presents numerical results comparing the accuracy, stability, and cost of the different algorithms.

Preliminaries on ellipses and ellipsoids
This section sets the stage for the comparison between different contact detection algorithms by collecting often recurring notation and definitions. By beginning with a compact but coherent introduction to the terms and expressions, we hope to make the similarities between the different algorithms quickly transparent.

Representation of ellipses
An ellipse E is the set of roots x = [x, y] T ∈ R 2 of a quadratic polynomial of the form where Q is a symmetric positive-definite (SPD) matrix in R 2×2 and c = [c x , c y ] T ∈ R 2 is the center of the ellipse. Formally written, an ellipse is defined as The coordinates x in which an ellipse is initially given will be referred to as the global coordinates, but there exists an isometry to a system of coordinates in which the geometry of E is especially straightforward. Indeed, a fundamental result of linear algebra states that for each SPD matrix Q, there exists an orthogonal matrix R, that is satisfying R −1 = R T and therefore in the form such that the matrix is diagonal with strictly positive entries, i.e., the eigenvalues of Q. An example of an ellipse is shown in Fig. 1 under the assumption that a ≥ b. The axes corresponding to a and b are called the semimajor axis and the semiminor axis, respectively. Accordingly, a and b are usually referred to as the semi-axes of the ellipse. Further including a translation to send the center c to the origin, we can introduce the local coordinates ξ = [ξ, η] T , with respect to which the ellipse consists in the set of roots of which can be recast in the classical form In this paper, f will be called the global geometric potential, or simply potential, of the ellipse, while f will be called the local potential. Clearly, the potential will always be a convex function with a minimum at the center c of the ellipse.
This definition allows one to interpret an ellipse E as the "circle" satisfying x − c 2 Q = 1 in the global coordinates. Eventually, when we consider the contact problem for two ellipses E i and E j , we may replace the subscript Q by the index i of the ellipse E i . Throughout this paper, the norm · written without a subscript will denote the usual Euclidean norm.
We now proceed with an explicit description of an ellipse as defined by (1). Let the matrix Q be given by where positive-definiteness is ensured by the conditions A > 0, B > 0, and det Q = AB − C 2 > 0. Then, in the global coordinates, the potential (1) is given as For convenience, we provide below the explicit relationships between D and R, and Q, namely An alternative form in which to express the potential is based on separating out the quadratic, linear, and constant terms.
Starting from (1), we find that the points x on an ellipse satisfy Introducing the potential f can then be rewritten in augmented matrix form as Another useful description is based on the parameterization of the ellipse in terms of a parameter t ∈ [−π, π[ such that all points given by lie on the ellipse. Using the mapping (4), the ellipse in the global coordinate system consists then of the points Certain algorithms for contact detection between two ellipses require the outward unit normal vector at a point x or, equivalently, at a point ξ , on an ellipse. From the definitions of the potentials f and f , see Eqs. (1) and (5), respectively, a simple calculation shows that the normal is given in global coordinates as or in local coordinates as Without delving into the explicit calculations, which can be found in several references [4,5], we note that the minimum radius of curvature is given by There exist several alternative descriptions of ellipses. For the sake of completeness, we mention here some of the most important descriptions even if they are not used in this review. The first one will nevertheless motivate an algorithm for finding initial guess points, to be detailed in Sect. 4.3. We recall that the focal points of an ellipse, f 1 and f 2 , are located on its semimajor axis, at equal distance from the center, and are explicitly given in local coordinates as The ellipse can then be defined as the set of points x that satisfy Moreover, it is possible to show that the normal vector at x generates a line that bisects the angle f 1 x f 2 . There is also a well-known description of an ellipse in terms of its eccentricity e = √ a 2 − b 2 /a ∈ [0, 1], with e = 0 corresponding to a circle [18]. Ellipses can be geometrically obtained as the cross section of the intersection of an inclined plane with a conic section, but this description of a 2-D ellipse requires three dimensions, thus making it less practical. Ellipses can also be described using mechanical means, such as the Trammel of Archimedes, the Tusi couple, or the ellipsograph of Benjamin Bramer [2,11,18]. The Steiner method for the construction of an ellipse is quite elegant but requires a discretization and is therefore not relevant to the continuous contact detection problem. In summary, ellipses possess a wealth of fascinating and unexpected properties, but none of these seem to be as useful as (1) or (5) when one needs to numerically estimate contact points.

Representation of ellipsoids
Similarly to ellipses, an ellipsoid E ⊂ R 3 is the set of roots to the potential: where Q is a 3 × 3 SPD matrix and c ∈ R 3 is the center of the ellipsoid. As in 2-D, there exists an orthogonal change of variable R ∈ R 3×3 , R −1 = R T , which diagonalizes Q such that D = R T QR. The eigenvalues of Q, i.e., the entries of D, are all strictly positive and are denoted by a −2 , b −2 , and c −2 , where the positive constants a, b, and c are assumed to be ordered as c ≤ b ≤ a. Using the change of variable (4), with ξ = [ξ, η, ζ ] T , one can write the local potential in its so-called local coordinate system (O, ξ, η, ζ ) as The positive constants a, b, and c are called the semi-axes of the ellipsoid.

Remark 1
Unlike in 2-D, the explicit form of the rotation matrix R can be obtained in several manners. We first note that an arbitrary ellipsoid is defined in terms of nine parameters: the coordinates of its center c = [c x , c y , c z ] T and the six entries of the symmetric matrix Q. However, since the matrix Q can also be written as RDR T , these six entries can also be identified with the three positive eigenvalues appearing as the diagonal elements of D, and the three parameters needed to describe a general orthogonal transformation R.
In other words, one needs to introduce three angles, each corresponding to an elementary rotation, and write R as the composition of the three rotation matrices in order to map the local coordinate system into the principal axes of the ellipsoid in the global coordinate system. The choice of the three angles is actually not unique and depends on the representation considered, for example the Euler rotations [19] or the quaternion rotations [21]. We will not describe these methods here and will simply assume that R, if necessary, is provided using one of the methods as An ellipsoid in the local coordinate system can be represented in terms of the spherical coordinates Using the mapping x = Rξ + c, the ellipsoid in the global coordinates system is therefore parameterized as As in 2-D, the outward unit normal vector at a point x on an ellipsoid is given in global coordinates by or in local coordinates by We conclude by observing that the gradient in 3-D is given by the same formula as (18), while the minimum radius of curvature is, assuming a ≥ b ≥ c,

Family of concentric similar ellipses and ellipsoids
Let E be an arbitrary ellipse or ellipsoid with potential Then, the family of concentric similar ellipses (d = 2) and ellipsoids (d = 3) associated with E consists of the sets (28) or as the roots of f r (x) We note that two ellipses or two ellipsoids within the same family, i.e., E(r 1 ) and E(r 2 ) with r 1 = r 2 , form a homoeoid, that is, the bounded region between E(r 1 ) and E(r 2 ). Moreover, E(0) reduces to the singleton {c}, while E(1) corresponds to E. Furthermore, for every point x ∈ R d , there exists a unique r ≥ 0 such that x ∈ E(r ) and the gradient of f r at x ∈ E(r ) is given by In other words, the outward unit normal vector to the ellipse/ellipsoid E(r ) associated with E at an arbitrary point x ∈ R d \{c} is given by We now provide some properties satisfied by the vector field n(x) associated with an ellipse E, which will be extensively used in Sect. 3. These properties will be expressed using complex multiplication and elements of projective geometry which we now recall. Let S 1 be the set of points of unit modulus in the complex plane C, which will be used to represent both the unit gradient field n and unit direction w. Given the points e iω and e iθ in S 1 , complex multiplication between the two points can be written as e iω e iθ = e i(ω+θ) , thereby representing the composition of two rotations.
In projective geometry, the real plane R 2 is imbedded into the compact space of all directions in R 3 using the association of [x, y] T ∈ R 3 with the direction is called the projective sphere SP 2 , not to be confused with the projective plane obtained after associating antipodal points in the projective sphere. The points at infinity are those corresponding to the directions [x : y : 0], thus forming a circle. In fact, given an affine map g : R 2 → R 2 , say for b ∈ R 2 and T a 2 × 2 matrix, then along the segment in direction w ∈ S 1 we can define a limiting direction This association is independent of c and parameterization t, thus leading to a well-defined map g ∞ : S 1 → S 1 according to This map is the restriction at infinity of the extension of g from the projective sphere to itself.

Lemma 1
Consider an ellipse E ⊂ R 2 centered at c ∈ R 2 whose unit vectors associated with the semimajor and semiminor axes are ξ and η, respectively, oriented counterclockwise. The vector field n given by (29) satisfies the following properties: (i) The vector field n(x) is well defined ∀x ∈ R 2 \{c}.
(iii) Given w ∈ S 1 , n(c + tw) is constant ∀t ∈ R + . (iv) The map is well defined and satisfies the following properties.
is uniform with respect to w ∈ S 1 .
Proof From (29), the vector field n(x) is the unit vector field associated with the gradient field Given that Q is SPD, ∇ f only vanishes at x = c. This proves Property (i). To demonstrate Property (ii), we observe that ξ and η are the eigenvectors of Q associated with the eigenvalues 1/a 2 and 1/b 2 , respectively; see (3). Hence, substituting x = c ± tξ into (33), we find which implies that n(x) = ±ξ . Similarly, substituting x = c ± tη into (33) shows that n(x) = ±η. Let w ∈ S 1 and consider the half-line supported by w, i.e., the set of points x = c + tw with t > 0. The gradients along the half-line are all positive multiples of the same vector Qw. Hence, the vector field n(x) is constant along the half-line, which proves Property (iii). This is illustrated in Fig. 2.
We now consider the proof of (iv) by first demonstrating (d). The bound (32) implies that the function is in fact the same as if we had taken lim n(x 0 + r w) and therefore does not depend on the origin x 0 of the segment x 0 + r w. For any x 0 = c and any direction w, For large and positive r , we have that As a matter of fact, this approximation can be made uniform in w for r sufficiently large. In other words, there exist constants R, δ, and C such that ∀r ≥ R and ∀x 0 ∈ R 2 satisfying Property (iv)-(a) follows immediately from ii), while Property (iv)-(b) follows from (iv)-(c). In particular, we observe that N (w) is the identity if and only if θ = 0, which, according to relation (31), occurs if and only if a/b = 1.

Contact detection problems
The purpose of this section is to introduce elementary notions of contact points and of their properties for pairs of ellipses. Lacking a common framework, much of the past work provided little indication of the connections between the different algorithms. This section is an attempt at filling this void by presenting a few precise definitions and results which will serve as a common thread in later comparisons of the different contact detection algorithms. The first element to comparing contact detection algorithms is to identify the mathematical contact detection problems which they are attempting to solve. In numerical analysis, this requires the demonstration of the existence of a unique solution to the contact detection problem, and often some stability bounds. This will be established for the three contact detection problems we have identified in the literature, but stability bounds and accuracy estimates for contact detection algorithms remain open.
Complete proofs will be provided in R 2 as some gaps in the theory remain in R 3 .
The Potential Contact Theory of Perram and Wertheim [40,41], based on the earlier work of Veillard-Baron [47], is formulated in algebraic terms that makes it one of the most mathematically complete theories to date. We will not be examining in detail the Perram-Wertheim theory because it does not provide estimates of separation/penetration distance [9], although it is widely used in static molecular dynamics where continuous force potentials are present. The continuous contact detection algorithm developed by Wang and his collaborators [6,51,52] is also well formulated mathematically, but it is too computationally expensive for the evaluation of the force in the quasi-static regime and therefore will not be described here. Nevertheless, we will make connections to those theories in Sects. 3.5 and 3.9. Unfortunately, the most computationally efficient algorithms are not always expressed with the same level of rigor, and the topic has suffered from some of this confusion.
In practice, every contact detection algorithm for ellipses should provide a single contact point, a single contact normal, and either a separation or a penetration distance. However, given two elliptical particles 2 E i and E j ⊂ R 2 , and two points judiciously constructed on each particle, say x i ∈ E i and x j ∈ E j , one could compute the distance between x i and x j as the separation or penetration distance and define the midpoint between x i and x j as the contact point (which should provide reasonable approximations of the contact properties in case of ellipses with small overlap). The contact normal could then be defined in terms of the segment joining x i to x j . For most of the algorithms, we shall describe in Sect. 5, this is precisely how the contact point and contact normal are actually computed.
Although the regularity of functions usually plays an important role in numerical analysis, all the objects studied below will be expressed as algebraic constraints between real analytic functions and hence will also be real analytic, which in practice means all functions will be infinitely differentiable away from a finite number of obvious singularities.
Estimating the contact point between all possible configurations of pairs of ellipses is in general not necessary in particle simulations and can lead to a number of degenerate cases, such as when the center of mass of one of the ellipses is inside the area of the second ellipse or when one ellipse is virtually penetrating completely through the other. We will therefore restrict ourselves to a limited number of configurations of the ellipses which we may encounter, illustrated in Fig. 3, thus avoiding degenerate contacts. Wang and his collaborators have produced a complete classification of all configurations, including degenerate contacts, using the theory of arrangements from computational geometry [27]. The first few definitions and lemmas below aim at characterizing such configurations, which we will refer to as near-perfect contact. When pairs of ellipses are in near-perfect contact, Theorem 7 will show that only a few cases need to be studied.

Intersection of ellipses
Their intersection set is defined as Before proceeding with an analysis of the intersection set, observe that it is, at least formally, computable as the solution to a minimization problem when I i j = ∅. Although not the only possible formulation, this section will show that its connection to other auxiliary problems makes it the richest formulation.

Lemma 2 Given two ellipses E i and E j such that
satisfyx i =x j :=x. According to this identification, the set of pairs (x,x) ∈ X i j correspond to the pointsx ∈ I i j .
Proof If the intersection is not empty, then there exists a point x ∈ E i ∩ E j for which the pair (x, x) ∈ E i × E j and the distance x − x vanishes. Any other pair of points (x i ,x j ) ∈ E i ×E j for which x i −x j reaches the minimum value, which we have shown is zero, must satisfy x i = x j . Thus, the single point x i = x j must belong to both E i and E j .
Intuitively, it is easy to imagine the different forms that the intersection set (see Fig. 3) may take, but it is less straightforward to give a complete and thorough description. Bézout's theorem [15,23] applied to the roots of two quadratic bivariate polynomials in R 2 states that the intersection set I i j can be either 1. empty: the ellipses are disjoint; 2. one point: the ellipses are in perfect contact (see Definition 4); 3. two, three, or four points: the ellipses intersect; 4. or an entire ellipse, if the two ellipses coincide.
We first state the following definition that will allow us to disregard trivial cases such as when one ellipse lies fully within a second ellipse or when two ellipses coincide. Although this condition will appear innocuous, it will play a key role in Lemma 11 and be a geometrical motivation for the numerical condition (60).

Definition 3 (Ellipses with non-penetrating centers of mass)
Two ellipses E i , E j ⊂ R 2 are said to have non-penetrating centers of mass (CoM) if the distances between the centers, evaluated in both the E i -and E j -norms (7), satisfy

Case of two disjoint ellipses
In this section, we consider the case of two disjoint ellipses E i and E j , i.e., I i j = ∅, that satisfy the non-penetrating CoM condition, see Definition 3. In this case, there is obviously no contact nor overlap, but one can estimate the distance between the two particles. Our objective in doing so is to find formulations of the separation distance that can be extended to the definition of distance, or contact point, when the ellipses are overlapping. The most obvious and straightforward formulation of the distance between two ellipses, which could naturally be applied to any pair of objects, is characterized in the following lemma.
The pair (x i , x j ) ∈ E i ×E j will be referred to as the minimum distance pair (MDP) of ellipses E i and E j .
Proof The existence of a unique pair (x i , x j ) ∈ E i × E j that minimizes distance x i − x j follows from elementary results in linear algebra [43], which we now present. Consider the convex and compact set E k := x ∈ R 2 ; f k (x) ≤ 0 , k = i, j, formed by E k and its interior. Since E i and E j are disjoint and have non-penetrating CoM, then E i ∩ E j = ∅. It implies that the set is also compact and convex. Hence, there exists a unique d ∈ D minimizing d , according to the Hilbert projection theorem. In fact, if d = x i − x j for (x i , x j ) ∈ E i × E j , then the pair must belong to E i × E j . If not, say x i ∈ E i \E i , then one could always find an ε, 0 < ε 1, such that Finally, we observe that the pair ( and thus, by virtue of the previous result, should lie on the boundaries of E i and E j . However, the points , form a straight segment and cannot lie on the boundary of E i (resp. E j ), since the ellipses are strictly convex. This shows that the pair Let f i and f j be the global potentials of E i and E j , respectively. In order to show that the unit normal vectors are opposite, we introduce the Lagrangian functional: where λ i ∈ R and λ j ∈ R denote the Lagrange multipliers associated with the constraints The solution (x i , x j ) to (42) is a stationary point of L and must satisfy or, equivalently, Combining those two equations leads to share the same or opposite direction. The identity (45) states that in the case that E i and E j are disjoint and satisfy the property of non-penetrating Similarly, the identity (46) shows that λ j < 0. Since both Lagrange multipliers are of the same sign, which allows one to conclude that (43) is satisfied.
It is worth noting that (43) represents a non-binding constraint as it is automatically verified by the solution to the minimization problem. Therefore, (42) can be recast as Remark 2 Although Lemma 3 is stated only for pairs of ellipses, the statement and the proof clearly also apply to pairs of ellipsoids.

Remark 3
In the case of two disjoint circles C i and C j , the solution pair (x i , x j ) to (42) is actually aligned with the centers of the circles, c i and c j . From this observation, one can reformulate the distance x i − x j as so that the minimization problem (42) can be recast as It follows that the minimization problem can be separated into the fully decoupled minimization problems In other words, the points x i and x j are the closest points on C i and C j to the centers c j and c i , respectively. Recalling that ellipses can be viewed as circles in their respective E-norms, one can actually introduce similar decoupled minimization problems in the case of ellipses.

Lemma 4 (Minimum Potential Pair) Let E i and E j be two disjoint ellipses with non-penetrating CoM with global potentials f i and f j , respectively. Then, there exists a unique pair of points
Moreover, following the convention (29), we have The unique pair (x i , x j ) ∈ E i × E j will be referred to as the minimum potential pair (MPP) with respect to the i-norm and the j-norm.
it follows that the two minimization problems in (48) are equivalent. The same reasoning implies that the two minimization problems in (49) are also equivalent. The demonstration of the existence and uniqueness to the minimization problems (48), or (49), is similar to the one given in Lemma 3. Consider the minimization problem where is compact and strictly convex. Then, it is well known that (52) has a unique solution, say x i . As we argued earlier, x i must in fact belong to the boundary E i and is unique because E i is strictly convex. The Lagrangian functional associated with the constrained minimization problem (48) is given by Since the solution x i to (48) is a stationary point of L i , it necessarily satisfies which, using the fact that the two ellipses are disjoint and have non-penetrating CoM, implies that the two normals at x i are in opposite direction. Indeed, as was explained in the proof of Lemma 3, it is easy to see that λ < 0, hence The relation (51) is shown in the same manner by introducing the Lagrangian functional L j associated with the minimization problem (49).
Since the relations (50) and (51) are satisfied at the points of the MPP (x i , x j ), they can each be added to the minimization problems (48) and (49), respectively, as non-binding constraints so that the two problems can be recast as and Remark 4 As we observed earlier after Lemma 3, the statement and proof of Lemma 4 apply equally well to pairs of disjoint ellipsoids.
Before proceeding with the other cases, we will make a few remarks on the solution to Problem (48). One classical approach for solving the constrained minimization problem proceeds by means of Lagrange multipliers, as seen earlier. However, the resulting problem could lead to several solutions as the nonlinear Lagrangian functional may have up to four critical points depending on the configuration and size of the ellipses. In other words, the solutions correspond to local minima and maxima of the potential function f j restricted to E i . This is exemplified in Fig. 4. In practice, all of the known methods identify all of the critical points and distinguish the global minimum by explicitly evaluating the distance at each critical point.
It is worth noting here that the non-binding constraint in the minimization problem (53) has the added benefit of yielding a Lagrangian functional with a unique critical point. Indeed, the constraint (50) is only satisfied at the global minimum. Alternatively, one may enforce the uniqueness of the critical point by considering the inequality constraint

Case of two ellipses in perfect contact
The case of perfect contact between ellipses with nonpenetrating CoM can be viewed as a limiting case of two disjoint ellipses. Therefore, we shall quickly verify that the previous results straightforwardly apply to this particular case.

Definition 4 (Perfect contact point)
Two ellipses E i , E j ⊂ R 2 are said to be in perfect contact if I i j consists of a single point and if they have non-penetrating CoM. That point is then called a perfect contact point.

Lemma 5 Let E i and E j be two ellipses in perfect contact
at point x c with non-penetrating CoM. Moreover, let n i (x c ) and n j (x c ) denote the outward normal unit vectors at x c to E i and E j , respectively. Then, the pair (x c , x c ) is the MDP and MPP of the two ellipses. Moreover, it holds that This shows that x c is the unique solution to (48), and, in a similar manner, x c is also the unique solution to (49). Finally, the relation (55) is clearly a consequence of (50) and (51)

Case of two ellipses with overlap
The goal of this section is to extend the definitions of MDP and MPP to pairs of ellipses with small overlap and to establish existence and uniqueness of those pairs. The description of the geometry at the intersection of ellipses with small overlap will therefore be a prerequisite for the analysis of a contact detection problem. Unfortunately, to make this intuitive notion precise we will need to introduce the co-gradient locus which formalizes the fundamental constraints (50)-(51) on the normals. In order to emphasize the key concepts and the structure of the theory, the proofs of the results will be delayed to Sects. 3.6 and 3.7. In contrast to the previous sections, most of the results below have only been established in R 2 and we will delay to Sect. 3.8 any speculation on extensions to R 3 .
The co-gradient locus is the set of points x where n i (x) ∝ n j (x), which includes points where (50)-(51) are satisfied. To the best of our knowledge, this set was first introduced in [45,46] where it was called the locus of common slope but was not characterized.

Definition 5 (Co-gradient function and co-gradient locus)
Given two ellipses E i , E j ⊂ R 2 , the co-gradient function is defined as The associated co-gradient locus is the set of all roots of the co-gradient function, i.e., In Sect. 3.6, we will provide a few alternative descriptions of h and its zero set H i j . The previous definition has an obvious analogue in R 3 . Building on Lemma 1, it is possible to prove the following elegant result illustrated in Fig. 5.
Theorem 6 (Co-gradient locus) Let E i , E j ⊂ R 2 be two ellipses with non-penetrating CoM. Then, the co-gradient locus H i j is a hyperbola and the two centers of the ellipses belong to the same branch of the hyperbola. The branch of H i j between c i and c j can be parameterized by a smooth injection Moreover, there exists a unique pair of parameters t k ∈]0, 1[, A long and detailed proof of this result can be found in Sect. 3.6, but there exists a shorter proof sketched out in Remark 6 of Sect. 4.2, based on mapping the pair of ellipses to a circle and an ellipse.
The curve H i j is an intrinsically defined invariant of the pair of ellipses, no matter whether they are in contact or not. The points γ i j (t i ) and γ i j (t j ) are therefore properties of only the pair and it would be natural to use γ i j (t i ) − γ i j (t j ) as a measure of closeness for pair of ellipses. Unfortunately, it is insufficient because it does not account for the relative size and orientations of the ellipses. The following numerical criterion will ensure a small overlap, in the sense of Theorem 7.

Definition 6 (Ellipses in near-perfect contact) Two ellipses
where This definition has geometric content that will become apparent in Sect. 3.6, particularly during the proof of the following theorem; see Sect. 3.7 for the proof.

and let t i and t j be defined as in Theorem 6. Then, the intersection
(1) if t i < t j , then I i j = ∅, i.e., the ellipses are disjoint; (2) if t i = t j , then I i j consists of a singleton, i.e., the ellipses are in perfect contact; (3) if t i > t j , then I i j consists of two distinct points, i.e., the ellipses have small overlap.
With this theorem, we have shown that the constraint (60) is sufficient to avoid degenerate contacts between ellipses, that is those with three or four contact points. The proof of this theorem also leads to an explicit description of the intersection; see Corollary 12. We are now in a position to introduce the extensions of MDP and MPP.

Theorem 8 Let E i and E j be two ellipses in near-perfect contact. Then, there exists a unique pair of points
In light of Lemmas 3 and 5, we can conclude that the correct formulation for the MDP, in all configurations of pairs of ellipses that are small perturbations of perfect contact, is the doubly constrained minimization problem (61). The proof of this fundamental result requires Theorem 6, Theorem 7, and Corollary 12; hence, it is deferred to Sect. 3.7. There is a similar result below for the MPP which shows that the general formulation of the MPP for all pairs of ellipses requires a condition on opposing normals, even though this condition was optional for disjoint pairs of ellipses.

Theorem 9
Let E i and E j be two ellipses in near-perfect contact with global potentials f i and f j , respectively. Then, there exists a unique pair of points (x i , x j ) ∈ E i × E j (known as MPP) satisfying the independent problems and This shows that x i = γ i j (t i ) must be the unique minimum of Problem (62). A similar argument applies to γ i j (t j ) ∈ E j ∩ H i j , thus completing the proof.

Contact potentials of Perram and Wertheim
The formulation of the MDP and the MPP has many similarities to the theory of contact potentials of Perram and Wertheim [41]. Although the purpose of contact potentials is only to identify contacts, and as such does not estimate the separation/penetration distance between particles, the mathematical similarities warrant a brief description. Perram and Wertheim developed functions, called contact potentials, whose sign could efficiently compute if two convex smooth particles were disjoint (> 0), overlapping (< 0), or tangent (= 0). This technique, extending earlier work of Vieillard-Baron [47], was widely adopted by statistical physicists studying phase transitions in specific models [9,40]. Today, simulations in statistical physics have also adopted even more efficient approximations to contact potentials, such as the Gay-Berne potential [16], for the generation of dense particle assemblies and fast computation of their thermodynamic properties. We note that contact potentials are equally well defined in both two and three space dimensions.
Consider two convex smooth particles, described by indexes i and j, whose interiors are defined by ζ i (x) < 0 and ζ j (x) < 0, respectively, and whose boundaries occur at ζ i (x) = 0 and ζ j (x) = 0, respectively. The global potentials of two ellipses satisfy these properties, but technically speaking the contact potentials do not require ζ i and ζ j to be of this form. The contact potential ζ i j for the two particles is the quantity and so it is clearly also defined as an optimization problem, just like the MDP and the MPP. For each fixed value of λ, the minimum with respect to z occurs at the critical point, i.e., the solution to This shows that for all λ ∈ [0, 1], the critical point z belongs to the co-gradient locus. Similarly, the maximum with respect to λ occurs at a value of z = z(λ) satisfying As Perram and Wertheim explained in their original paper, ζ i j corresponds to the value of at the point x ∈ H i j where equal scalings of the original surfaces, in the sense of (65), are in perfect contact, according to Definition 4. More details about the geometry of contact potentials can be found in an article of Paramonov and Yaliraki [39]. The development of the theory of contact potentials is largely expressed in terms of explicit algebraic constraints, and so most questions surrounding existence and uniqueness are rather easy to deal with. In fact, many of the key results were already demonstrated in an appendix to the original paper by Vieillard-Baron [47].

The co-gradient locus
The analysis of the previous sections has highlighted the key role that the co-gradient locus, see Definition 5, plays in characterizing intersections for ellipses in near-perfect contact, and even in solving for the MDP and the MPP. The purpose of this section is to make a thorough study of H i j in R 2 , including when the hyperbola degenerates. We begin by offering a few equivalent definitions of the co-gradient locus that will be useful later on.
In 2-D, the cross-product defining the co-gradient function is interpreted as a cross-product in 3-D between the gradients in the 2-D plane. It therefore results in a 3-D vector with only one nonzero component along the z-axis, which is given by the scalar function We shall consider this definition of the co-gradient function when dealing with ellipses rather than the vector-valued h in (56). Introducing the anti-symmetric matrix, the co-gradient function can also be written as where we have used the fact that Q i is symmetric. We immediately observe that h is a quadratic polynomial in x and that the centers c i and c j of the ellipses belong to H i j . If the product Q i AQ j was symmetric, then the determinant of the product could immediately tell us the geometry of the cogradient locus. Unfortunately, the detailed characterization of H i j presented in Theorem 6 requires significantly more work.
A second characterization can be made by normalizing the gradients in (56). Recalling the definition of the unit normal vectors (16), i.e., then the normalized co-gradient function iŝ The scalar component in the is the angle between n i and n j , well-defined modulo 2π . Mimicking the definition (68) of the z-component ofĥ, the co-gradient locus in 2-D is simply the set of roots of The angle η i j (x) can also be defined by identifying n i and n j with unitary complex numbers, so that, using complex multiplication The roots ofĥ correspond to η i j (x) = mπ , m ∈ Z. We note that at infinity, relation (72) can be rewritten as N j (w) = e iη i j (w) N i (w), where N i and N j are the maps of Lemma 1, that is the angle η i j = η i j (w) only depends on the direction w. Eventually, we will show that the angle η i j must belong to ] − π, π[, and hence is well defined. We now proceed with the proof of Theorem 6.
Proof The proof will show that the co-gradient function h, which is a quadratic function according to (68), possesses four roots at infinity. This will imply that the roots H i j ⊂ R 2 form a hyperbola because ellipses, parabolas, and hyperbolas possess, respectively, 0, 2, and 4 roots at infinity on the projective sphere. Afterward, we will argue that a single branch of the hyperbola must cross both centers of the ellipses, thereby justifying the existence of the parameterization. The majority of the analysis will be performed on a pair of ellipses in a generic configuration but this will require us to begin the proof with a lengthy study of different degenerate configurations. Bivariate quadratic polynomials have roots that can degenerate to either of the following configurations: two intersecting lines, two parallel lines, a line with a second line at infinity, two coincident lines, or a single point. The last option will never occur because h already vanishes at the centers c i = c j . The analysis below will show that H i j always contains at least four points at infinity and hence cannot be formed of two parallel lines or two coincident lines.
The first configuration we study assumes that c j belongs to the axis ξ i and that the principal axes of E j are aligned with those of E i , although the argument will also work if c j belongs to the axes η i and η i = ξ j . Under these conditions, for all t, s ∈ R, Property (ii) of Lemma 1 shows that Hence, every point of the axis c i +tξ i belongs to H i j . Furthermore, the fact that the axes are aligned and Property (iv)-(a) implies that or, in other words, h possesses roots at infinity in the direction w = ±η i = ±η j . In this configuration, the co-gradient locus has four roots at infinity and contains a line and hence is either two intersecting lines, or a line with a second line at infinity. If the two ellipses have the same aspect ratio, then relation (31) of Lemma 1 states that N i (w) = N j (w), for all w ∈ S 1 . In other words, H i j contains the line at infinity. On the other hand, when the aspect ratios are different, the same relation shows that N i = N j and thus the co-gradient locus must be formed of two intersecting lines. Fig. 6 Illustration of the proof of Theorem 6. The ellipses E i and E j are in a configuration with c i at the origin and E i is aligned with horizontal axis. The circle with radius r is large enough that all normals are external on the ellipses, i.e., n i · n j ≥ 0 Consider now the case where c j does not belong to either principal axis of E i but continue to assume that the principal axes of both ellipses are aligned, say ξ i = x j and η i = η j . Property (iv)-(a) of Lemma 1 tells us that and hence there are at least four roots at infinity. Again, if the ellipses have the same aspect ratio, then N i = N j and the co-gradient locus is degenerate and contains a line at infinity. Otherwise, the co-gradient locus is a hyperbola, which or may not be degenerate; see Corollary 10 for more on this issue. The general configuration on which we will focus the remainder of our attention assumes that the principal axes of the two ellipses are not aligned, whether or not either center belongs to the axis of its brethren. In this case, we observe that the straight line c j + tξ j , t ∈ R, for |t| sufficiently large, belongs to two opposing quadrants: either When the second case occurs, the line tξ i crosses the opposing quadrants [ξ j , η j ] and [−ξ j , −η j ]. Hence, the second case can be brought into the first configuration by translating E j to the origin and exchanging the indices i and j. Figure 6 illustrates this configuration, after assuming a translation and a rotation sending c i to the origin and the axes of E i over to the usual Cartesian axes.
The map (30) associates with each direction w ∈ S 1 , the unique normals N i (w), N j (w) on the line at infinity. We may then measure the angle η i j = η i j (w) between the normals at infinity using the relation (72), rewritten here using complex multiplication as In this last identity, positive or negative angles correspond, respectively, to a counterclockwise or clockwise rotation when rotating N i toward N j . It is essential to observe that the angle η i j is well defined within ] − π, π[ because Property (iv)-(c) of Lemma 1 states that the rotation θ from w to either N i (w) or N j (w) is strictly bounded |θ | < π/2, and hence the rotation from N i to N j must be by an angle η i j strictly less than π in absolute value. As w moves counterclockwise around S 1 starting at ξ i , the configuration we have chosen, as shown in Fig. 6, implies that we will encounter in order the directions We will focus on demonstrating that η i j possesses a root inside the arc [ξ j , η i ] ⊂ S 1 , but similar arguments will show that there are at least three other roots, one in each of the three arcs Each root of η i j corresponds to equal normals and hence to a root of h, thereby demonstrating that H i j is a hyperbola. Consider the principal axis for E j centered at c j = c i , i.e., take w = ξ j . Then, estimate (32) and Property (iv)-(d) of Lemma 1 show that the unit normal n j (c i + r ξ j ) converges to N j (ξ j ) = ξ j as the radius r increases. Furthermore, if ξ j = e iσ j ξ i for σ j ∈ [0, π/2[, then Property iv)-(c) states that Given that tan σ j > 0 and a i /b i > 1, we find that θ i > 0 and therefore must belong to [0, π/2[. Using these facts and the estimate (32), we find −π, π[ and changes sign as the direction w varies from ξ j to η i , there exists a direction w in ]ξ j , η i [ where η i j vanishes and both normals are equal. The same argument is applied to the three other intervals, which shows that there are four directions at infinity where the normals coincide, i.e., the co-gradient function h has four roots at infinity. To conclude the proof, we need to show that a single connected branch of the hyperbola crosses the centers of the two ellipses. It is easy to verify that c i and c j belong to H i j by substituting directly into (68). Consider the function α(x) = n i (x)· n j (x) for all x belonging to the smooth affine variety H i j . By construction, α only takes on the values ±1 and is ill-defined at both centers. Since α is continuous on H i j \{c i , c j }, it must have constant values over each connected component of H i j \{c i , c j }. Furthermore, for points x far from the centers, the normals are related by n j = e iη i j n i where η i j ∈ ] − π, π[ and therefore the dot product must take positive values, i.e., α(x) = +1. If the points c i and c j belong to different branches of the hyperbola H i j , then each connected component of H i j \{c i , c j } reaches infinity and α must be equal to +1 everywhere. Yet, in the neighborhood of a center, say c i , the normal n j varies smoothly, and along the tangent to H i j at c i , property (ii) of Lemma 1 states that the normals n i are equal and opposite on both sides of c i . Hence, the function α must take opposite values, i.e., +1 and −1, when passing through the center c i along the branch of the hyperbola. This would contradict the conclusion that α is identically +1 on H i j \{c i , c j }, following from the hypothesis that c i and c j belong to different branches of H i j . The existence of the parameterization γ i j is a trivial consequence of the fact that both centers belong to the same branch of the hyperbola. The fact that H i j is a hyperbola implies that the portion γ i j ([0, 1]) ⊂ H i j can only cross each ellipse E k , k = i, j, once when connecting both centers, thereby demonstrating the uniqueness of t k , k = i, j.
The next result is included for the sake of completeness because it clearly shows that properties of the pair of ellipses are reflected in the geometry of H i j . In this sense, it indicates that additional study of the co-gradient locus might prove worthwhile, particularly the study of the second branch. Nevertheless, this result will only be used in the construction of an example later in this section. Proof The analysis in Lemma 6 has already shown that the degenerate quadratic cannot be formed of a single point (c i , c j ∈ H i j ), two parallel lines or two coincident lines (H i j at infinity always has ≥ 4 points). The only two remaining possibilities are those mentioned in the statement of the corollary. To complete the proof, we will first deduce the conditions required for H i j to contain a line at infinity. Afterward, we will identify the conditions under which the co-gradient locus contains at least one line. If the ellipses have aligned axes, then relation (31) clearly implies that the co-gradient locus possesses a line at infinity if and only if the aspect ratios are the same. Suppose now that the axes are not aligned but that the co-gradient locus still possesses a line at infinity. We will show that these hypotheses lead to a contradiction. Considering the axis directions as elements of S 1 , suppose that ξ j ∈ ]ξ i , η i ], although similar arguments will apply if it belongs to ]η i , −ξ i [. Given that H i j possesses a line at infinity, we have that N i = N j and in particular

Corollary 10 Two ellipses
The map N i cannot be the identity map, or else both ellipses would be circles and hence have aligned axes. Hence, for The first relation in (74) for ξ j = e iσ ξ i gives us that is θ = 0. The implicit relationship (31) between θ and σ shows that θ can vanish only if a i /b i = 1. Repeating the same argument with the second relation in (74) proves that a j /b j = 1. In conclusion, if the axes are not aligned but H i j possesses a line at infinity, then both ellipses are circles. We now attempt to determine the most general conditions under which the co-gradient locus can degenerate to a pair of intersecting lines. If we examine the values of the continuous function α(x) = n i (x) · n j (x) along H i j , we notice that it only take on the values ±1. During the proof of Theorem 6, we observed that α changes sign as we cross either center but that α on co-gradient locus always took the value α(x) = 1 when x was sufficiently large. This implies that if H i j degenerates to two intersecting lines, then the two centers cannot belong to different lines, and when they do, the second line cannot cross the first between the two centers.
Consider the line connecting both centers, which can be parameterized as In the previous identities, the angles were functions θ k = θ k (σ k ) for w = e iσ k ξ k , k = i, j, according to the implicit relations Since N i (w) = −N j (−w), we conclude that Each ellipse is parameterized in a space of dimension 5, two for each center c k and three for each matrix Q k . These 10 parameters determine a k , b k , θ k,0 , σ k , and θ k ; hence, the three identities (76) and (75) determine a subspace of codimension 3 where H i j is degenerate. In other words, Sard's theorem [33] states that H i j degenerates only over a subset of measure zero in the space of parameters for E i and E j .

Remark 5
We will present a degenerate co-gradient locus formed of two intersecting lines but for a pair of ellipses whose principal axes are not aligned. This example is instructive in that it goes beyond the degenerate examples identified during the proof of Theorem 6 and demonstrates that the conditions (75) and (76) are non-empty. Furthermore, this example is new to the literature and could be used for code verification.
Consider the ellipse E i described by its geometric potential The point x i = [ √ 3/2, 1/ √ 2] T belongs to E i and the line passing through the origin and x i forms an angle of π/6 with the horizontal axis because tan(π/6) = 1/ √ 3.
Using Formula (16) shows that the gradient to E i at x i is given by A simple calculation shows that the angle measured counterclockwise between x i and ∇ f i (x i ) is also π/6. This conforms with relation (31), Our objective is to construct an ellipse E j whose center c j belongs to the line and whose normal N j in the direction −x i is opposite to the normal N i in the direction x i . To keep this as simple as possible, we will describe the ellipse in its local coordinates and we introduce the point on the ellipse. With respect to the ξ j axis and measured counterclockwise, the line through the origin and x j forms an angle of 5π/4. The gradient to E j gradient at x j is and it is easy to see that it forms an angle of π/6 with x j . Hence, relation (31) is again satisfied using the fact that tan(5π/12) = 2 + √ 3. The ellipse E j can be translated so that c j is located along the line , say in the first quadrant, and then it can be rotated by −π/12 = π/6 − π/4 so that the point x j also falls on the line. Afterward, the axis of the two ellipses will no longer be aligned, but their normals will be opposite everywhere on the line between the two centers. The line will therefore be a part of the co-gradient locus. We conclude by remarking that this example clearly satisfies the three conditions for degeneracy of Corollary 10, that is formulas (75) and (76).

Characterizing near-perfect contacts
This section will demonstrate the existence and uniqueness of the MDP for pairs of ellipses in near-perfect contact, see Theorem 8, after proving the characterization of near-perfect contact in Theorem 7. The proof of Theorem 7 includes an explicit characterization of the intersection, summarized in Corollary 12, which is very useful when studying the MDP. The results in this section are again stated in R 2 , but the extensions to R 3 are rather straightforward.
Rather than attempting directly to study the intersection of two ellipses, we begin with a first-order approximation with two circles. Once equipped with explicit estimates with pairs of overlapping circles, we will eventually argue that the intersection of pairs of ellipses can be approximated by intersections of pairs of circles. Although the detailed results of this section were not all necessary for the study of the MPP, as the proof of Theorem 9 shows, they are necessary for MDP. In fact, we expect that the tools described in this section will find other uses.

Lemma 11
Consider a pair of overlapping circles C i and C j with non-penetrating CoM. If D i and D j are the closed discs bounded by the respective circles C i and C j , then the intersection is the union of two closed domains A + i j and A − i j , Proof We begin by simplifying the geometry through a translation and a rotation sending the center c i to the origin, and the second center c j to (c j , 0) along the positive x-axis. Corollary 10 states that the co-gradient locus is the x-axis with a second branch at infinity. The potential for the circles is and the normalized distance to the centers is We now proceed to detail the intersection between the two discs. The circles are overlapping, hence c j < r i + r j . The condition of non-penetrating CoM implies that max{r i , r j } < c j . The intersection can be characterized as This allows us to define the triangular domain Our objective is now to show that the maps are bijective diffeomorphisms. Given a point (λ i , λ j ) ∈ T , we will find a point x ∈ A + i j , or A − i j , such that For x = (x, y), the map is If we isolate the x coordinate, we find Substituting back into the equation for λ 2 i and simplifying This equation possesses two solutions for y, corresponding to either the map ϕ + or ϕ − . Taking λ i = λ j = 1, we can quickly verify that there are only two solutions at the intersection of the two circles. From the construction, it is clear that ( The previous calculation shows that the maps ϕ ± are bijective between A ± i j and T . To show that they are diffeomorphisms, we compute the differential by first observing that If one recalls the definition (68) of h, one can conclude that In the case of two circles, a simple calculation shows that which only vanishes along the co-gradient locus, that is, along part of the boundary of A ± i j . The proof of the next theorem will show that the condition of near perfect contact implies that the pair of ellipses satisfy two geometric properties. To explain these conditions, we introduce D + k (r ) the disc of radius r tangent to E k at γ i j (t k ) and whose interior is disjoint from Similarly, let D − k (r ) be the disc of radius r tangent to E k whose interior lies inside E k . Intuitively, D + k (r ) (D − k (r )) is the disc tangent to E k placed on the outside (inside) of E k , see Fig. 7a. We will also be using ρ k and ρ k the minimum and the maximum of the radius of curvature of the ellipse E k , respectively.
Definition 6 implies that (1) the discs D − i (ρ i ) and D − j (ρ j ) have non-penetrating CoM; and (2) the portion of H i j inside the intersection of the two ellipses is entirely inside the intersection Both conditions are generalizations of what was observed in the proof of Lemma 11. We now proceed with the proof of Theorem 7.
Proof This will not be a complete proof of Theorem 7 but an attempt at explaining how the condition is related to the two geometrical properties in the previous remark, as well as the maps ϕ ± describing the intersection. We begin by studying the neighborhood of the point γ i j (t i ) ∈ E i . The tangent to the curve γ i j at γ i j (t i ) is not necessarily in the same direction as the normal n i (γ i j (t i )), but certainly not perpendicular, and we note the defect as We now establish a condition under which the tangent line . Consider a translation of γ i j (t i ) to the origin, then apply a rotation to send the tangent line to E i at t i to the horizontal x-axis, as in Fig. 7b. We observe that we can reduce the analysis to showing that a curve crossing the origin at an angle η with respect to the vertical axis remains inside the discs tangent to the x-axis centered at the origin.
Recall Formula (18) for the smallest radius of curvature ρ i := b 2 i /a i , and the parameterization of the points on the boundary of D − i (ρ i ), centered at (0, −ρ i ) below the horizontal axis: The smooth curve γ i j remains close to its tangent line sγ i j (t i ), ∀s ∈ R, as long as t stays close to t i . We now compute the length of the portion of the tangent line inside D − i (ρ i ), which by symmetry will be of the same length inside D + i (ρ i ). The tangent line forms an angle η with the vertical axis. Let p = 0 be the unique point on D − i (ρ i ) where the tangent line crosses. We observe that an equilateral triangle is formed between the origin 0, the center (0, −ρ i ), and the point p with two sides of length ρ i and two angles of measure η. The length we are looking for is the length of the side of this equilateral triangle opposite to the center (0, −ρ i ). Based on this geometry, it is easy to verify that the length of the portion of the tangent line inside The curve γ i j (t) will remain close to its tangent line as long as the parameter t stays close to t i , with the deviation being proportional to |t −t i | 2 and the curvature of γ i j . Requiring that the parameterized distance |t i − t j | along the curve (a) (b) Fig. 7 (a) Illustration of disks D ± k (ρ k ) at points γ i j (t k ) with k = i, j for a pair of ellipses E i and E j . (b) Transformation of the two disks D ± i (ρ i ) with γ i j (t i ) sent to the origin and the tangent line to E i at t i aligned with the horizontal x-axis γ i j be small when compared to the length above implies that the curve cannot exit D + i ∩ D − i and hence the portion of H i j with t close to both t i and t j crosses E i only once. The minimum over i and j in (60) also constrains the co-gradient locus, at least the part inside D − j (ρ j ) ∪ D + j (ρ j ), to contain only a single point from E j . It is also clear that by taking t j − t i sufficiently small and negative, then the two discs D − i (ρ i ) and D − j (ρ j ) will satisfy the non-penetrating CoM condition.
We can now distinguish the following three cases: • If t i = t j , then the definition of the co-gradient locus implies that both ellipses have opposing normals and that the two ellipses are in perfect contact. • Suppose now that both t i < t j and (60) are satisfied.
The portion of the co-gradient locus between γ i j (t i ) and , which is outside of both ellipses. There exists a unique ellipse E j (r j ) which is tangent to E i at γ i j (t i ) and because γ i j (t i ) is outside of E j , r j must be greater than 1. Since both E i and E j (r j ) are convex and tangent, E i is disjoint from E j (r j ) and the set E j is strictly bounded by E j (r j ). These remarks imply that the ellipses are disjoint when t i < t j .
• Suppose now that both t j < t i and (60) hold. The objective is to show that E i ∩ E j contains exactly two points, which is fundamentally no longer a question only about a neighborhood of γ i j (t) and t ∈ [t i , t j ] but a global question about a region. The question will be answered by providing a detailed description of the region at the intersection E i ∩ E j . The characterization of the intersection E i ∩ E j will be a simple corollary of the description of E i ∩ E j . We shall show that there exists a diffeomorphism of a triangle toward each half of the domain E i ∩ E j \H i j .
We begin by identifying a subdivision of E i ∩ E j . Following the definition of h in (66), we can construct and deduce that Intuitively, the normalized co-gradient function (71) can be used to uniquely define the angle η i j (x) between the two normals in a neighborhood of H i j according tô and the convention that we measure increasing η i j when rotating counterclockwise from n i to n j . Hence, for all x in a neighborhood of H i j = A + i j ∪ A − i j , there is a well-defined angle η i j (x) with values near π . Since the scalar function h is the z-component of the cross-product n i × n j , the two regions A + i j and A − i j correspond to the regions where η i j (x) is, respectively, less than π and greater than π , see  As introduced earlier in Sect. 2.3, each point x ∈ R 2 belongs to unique scaled ellipses E i (λ i (x)) and E j (λ j (x)) where the real-valued functions are smooth away from the centers c k . The gradients of these functions are multiples of the gradient because The values of λ i and λ j along the co-gradient locus belong to the intervals [r i , 1] and [r j , 1], respectively, where r i := λ i (γ i j (t j )) and r j := λ j (γ i j (t i )). More precisely, along the segment γ i j ([t i , t j ]), the parameter λ i (γ i j (t)) is increasing and λ j (γ i j (t)) is decreasing and therefore we can parameterize λ j as a function of λ i along the segment, say λ j (λ i ). From this, we can define a triangular domain and two smooth maps The determinants of the differentials of these maps are hence the maps ϕ ± are of rank 2 except along H i j . This implies that each map is locally bijective, at least in the interior of A ± i j . The analysis presented so far shows that ϕ ± is a diffeomorphism only in a neighborhood of H i j . A complete proof would require a proof that λ j be monotone along the boundary of E i in A ± i j , and vice versa for λ i . This would demonstrate that ϕ ± is bijective along the boundaries. Given that D − k (ρ k ) is tangent to E k , it is intuitively clear that if Lemma 11 holds for both discs, then the maps ϕ ± in the case of ellipses should also characterize the intersection E i ∩ E j .
The proof of Theorem 7 given above also provides the following explicit description of the intersection.

Corollary 12 Consider two ellipses E i , E j ⊂ R 2 in nearperfect contact. If E i and E j are the closed regions bounded by the respective ellipses, then the intersection is the union of two closed domains
The proof of existence and uniqueness of the MDP for pairs of ellipses in near-perfect contact, namely Theorem 8, requires both Theorem 7 and Corollary 12. Although the result is intuitively convincing, the proof is long and will not be included in this review because the formulation of the MDP and the structure of the theory is, in the authors opinion, more important than the details of the individual proofs.

Extension to ellipsoids
The objective of this paper is to review and compare algorithms for the rapid and accurate estimation of separation/penetration distance for pairs of ellipses and ellipsoids in the quasi-static regime. The main mathematical results are Lemma 1, Theorem 6, and Theorem 7, which were all stated and demonstrated in 2-D, so it is natural to enquire as to what can be said in R 3 . In this section, we will briefly suggest how these three results could be extended to 3-D. Before proceeding, we remark that Definitions 4 and 5 have obvious extensions, that Lemmas 3 and 4, as well as Theorem 5, have identical statements and proofs in 3-D, and that in practice no change is required to the algorithms in Sect. 5.
The most difficult result to extend to 3-D is Theorem 6, and it should be studied with the techniques of algebraic geometry [20,23]. Lemma 1 was only a preliminary result for the proof of Theorem 6 and provided tools to circumvent real projective geometry; hence, we will discuss only Theorem 6. As far as Theorem 7 is concerned, the extension to 3-D is straightforward assuming that H i j has been characterized. In fact, our initial proof of Theorem 7 was carried out in 3-D and the key map ϕ is slightly easier to study because the intersection E i ∩ E j is path-connected (but not simply connected). In R 3 , the additional parameter in the intersection E i ∩ E j parameterizes a circle S 1 whose tangent is given by The statement in 3-D of Theorem 6 concerns the cogradient locus H i j , which is now the set of common roots of each of the three components of the cross-product (56), each of which can be written roughly in a quadratic form similar to (68). Each component defines a two-dimensional subvariety that intersects the sphere at infinity in real projective space RP 3 along a curve. We conjecture that, as we did in 2-D, we can characterize the common zeros in R 3 by identifying the isolated intersection points of the three 1-D curves at infinity. Each ellipse will define three planes corresponding to each pair of its orthogonal axis, where the normals take on known values. These six curves on the sphere at infinity will provide a triangular subdivision of the sphere and the existence of a root to (56) in each triangle will be determined by examining the signs of the components of h along the edges and nodes of the triangle.

Relationship to time-dependent contact detection
For pairs of rapidly moving and/or rotating ellipses/ellipsoids, estimating the contact point and the penetration distance requires anticipating future positions. More specifically, the future contact point might be very different from the MDP at a given time when the relative velocity of the two particles is large and their surfaces are close. In this section, we show that the constraint of the co-gradient locus H i j appearing in the MPP can be interpreted as enforcing displacements of the ellipses along a specific trajectory. In other words, we can interpret H i j as a specific particle trajectory bringing the two ellipses in contact. As we mentioned in the Introduction, the purpose of this review is not to study time-dependent collision detection for which there exist particular techniques, such as the continuous contact detection method of Wang et al. [6], but to highlight similarities that might allow one to formulate analogous results to those in Theorems 8 and 9 for time-dependent problems. Consider two ellipses E i , E j that are disjoint but in nearperfect contact. Consider the parameterization γ i j of the curve H i j going from c i = γ i j (0) to c j = γ i j (1), and fix t ∈ ]t i , t j [ corresponding to γ i j (t) in the exterior of both ellipses. Along the directed line c j γ i j (t), there is one unique intersection point p j (t) ∈ E j , which depends analytically on t, and the translation is defined by Remark that γ i j (t) − c j and γ i j (t) − p j (t) are colinear; hence, the displaced ellipse T t E j has the same normal along the line T t c j to T t p j (t) as it had along the line connecting c j to γ i j (t). This implies that T t will map p j (t i ) to γ j (t i ) and that at γ i j (t i ), the normal to T t E j will be the same as n j for E j . The definition of H i j then implies that the image of T t i E j will be in perfect contact with E i .
We have therefore shown that as t goes from t j to t i , the transformation T t brings E j in contact with E i at γ i j (t). This transformation does not correspond to a free fall trajectory for E j , because on the one hand, each translation T t preserves the orientation and therefore the physical displacement could not have applied torque. Yet, the displacement of the center of mass T t c j is determined (indirectly) by the points on a hyperbola in any rotated frame of reference, and therefore does not follow a parabola. Hence, T t does not describe a physical displacement. Nevertheless, this construction provides an additional interpretation of H i j as a canonical trajectory associated with each pair of ellipses.
Finally, we would like to highlight the research of Wang and his collaborators done for rapidly moving particles. That research did not attempt to measure separation/penetration distance but only to identify whether the ellipsoids are disjoint or overlapping (or in perfect contact). In 2001, Wang, Wang and Kim [52] found an algebraic condition for the separation of two ellipsoids. They demonstrated that a certain characteristic polynomial p i j of degree 4 had two positive real roots precisely when the ellipsoids were disjoint. In [51], the authors showed that the earlier algebraic condition could be used to identify a plan separating the two ellipsoids. The use of a separating plane was shown to be criteria that could be reused for several consecutive timesteps, without additional computation, if the linear and angular velocities of the particles were not too large, and hence lead to significant savings. The work in [26] proposed to apply an improved criterion [3] for the determination of the sign of the roots of the characteristic polynomial p i j from the original work [52]. In the spirit of Descartes' rule of signs, this extension involved using determinants of submatrices of the ith Sylvester-Habicht matrix, which is a family of reduced resultant matrices. This technique produced increased computational efficiency for the contact detection problem, but still without a computation of the distance between the ellipsoids.

Initializing the positions of two ellipses for contact detection
Some algorithms to be presented in Sect. 5 share two common features. The first is that the problem can be simplified by introducing normalized coordinates where one of the ellipses has become a circle. The second common feature is that they may require an initial guess of the contact point, and obviously, the efficiency of an iterative algorithm will be highly dependent on the quality of that initial guess. In this section, we review two mapping steps; see Sects. 4.1 and 4.2. The mapping will simplify the later presentation of the algorithms, and one of our proposed techniques for guessing the contact point (see Sect. 4.3) is done in one of the normalized coordinate systems we introduce below. We remark that the estimation of the contact point is rarely addressed in the literature and, to the best of our knowledge, the focal point estimate was first presented by the authors in [29].

Mapping of (E i , E j ) into a unit circle C i centered at origin and an ellipse E j
The first mapping was suggested by Ting et al. in [46]. Let E i and E j be two arbitrary ellipses defined by where the diagonal matrices D i and D j and the rotation matrices R i and R j are as defined in (3) and (2), respectively. The mapping consists in transforming one of the two ellipses, say E i , into the unit circleĈ i centered at the ori-gin and the other ellipse E j into the ellipseÊ j ; see Fig. 9. The mapping that transformsĈ i into E i consists of the transformation D −1/2 i , followed by the rotation R i of angle θ i and the translation c i . We thus have Indeed, it is straightforward to check that the equation of the circleĈ i in the coordinate system (Ô,x,ŷ), using (82) in (80), thus readŝ The equation of the new ellipseÊ j , following the mapping of ellipse E j , is obtained by substituting (82) in (81). We find thatÊ j in the (Ô,x,ŷ) coordinates are the roots of where we have introducedĉ j satisfying The equation forf j reduces tô wherê For the sake of completeness, we remark that the coefficients ofQ j , i.e., can be found explicitly aŝ and the parameters {â j ,b j ,θ j } associated with ellipseÊ j , so thatQ j =R jD jR T j , can be recovered by identification from Formula (10). Fig. 9 Mapping of (E i , E j ) into a unit circleĈ i centered at the origin and an ellipseÊ j The following lemma explains that the minimization problem in the new coordinates is related to a minimization in the original coordinates.

Lemma 13 In the coordinate system (Ô,x,ŷ), the Euclidean distance is the same as the distance with respect to the E inorm in the original coordinates (O, x, y).
Proof For every x, its imagex is computed as We then have, for arbitraryx andŷ,

Mapping of (E i , E j ) into a unit circleC i and an ellipseĒ j centered at origin
Džiugys and Peters [12] suggested another mapping such that the ellipseÊ j introduced above is now positioned in its local reference system denoted here as (Ō,x,ȳ). We will describe this new mapping as the previous mapping followed by a rotation and a translation sending the previous ellipsê E j to the origin with axes aligned with (Ō,x,ȳ). The circlê C i , previously at the origin under the mapping of Sect. 4.1, is now shifted around the ellipse centered at the origin, see Fig. 10.
The mapping amounts to considering the transformation that mapsx intox in R 2 by the rotationR j and translation c j : so that circleĈ i , located at the origin, is now mapped into circleC i with centerc i as The equations of the two ellipses in the new coordinate system (Ō,x,ȳ) are then given bȳ In other words, the original ellipse E j is now transformed into the ellipseĒ j centered at the origin and the ellipse E i is mapped into the unit circleC i centered atc i using the global transformation obtained by combining (82) and (86) where we have used (84) to obtain the last expression.

Remark 6
The transformation of an arbitrary pair of ellipses into an ellipse in its local coordinate system and a circle can be useful to simplify the mathematical analysis of the pair of ellipses. As an example, we can easily show that the cogradient locus in the given configuration (87) is a hyperbola.
Recalling the co-gradient function (68), the equation of the co-gradient locus reads in this case Developing the above equation leads to This is actually the equation of a hyperbola in the case that a i = b i . Indeed, using classical formulas, the center of the hyperbola is given by and, using the change of variables ξ = x −x h and η = y− y h , the equation can be recast as which is the equation of an hyperbola. In the case that a j = b j , the locus reduces to a straight line passing through the origin (i.e., the center of ellipse E j , which is a circle here) and the center of circle C i .
In the next section, we shall see that the mapping allows one to define accurate initial guess points for the solution of the minimum potential algorithm (GPA).

Initialization of contact point algorithms
The purpose of this brief section is to discuss possible initial guesses for contact points, taking for granted that the various algorithms will often rely on iterative methods to produce the final estimate of the contact point. In practice, the accuracy and efficiency of iterative methods are closely tied to the quality of the initial guess. In fact, after each timestep of the DEM it is common to use the estimate of the contact point at the previous timestep as an initialization at the current timestep, and so the initial guess may only be required when two ellipsoids are first identified as a pair with a potential contact during the broad phase search. Nevertheless, this question of fundamental importance, that is rarely addressed, should at least be briefly discussed in a thorough review.
Let us consider an ellipse E in its local coordinate system and a unit circle C in near-perfect contact. (We dropped indices i and j here for the sake of simplification in the notation.) The goal here is to identify a reasonable approximation p of the contact point x ∈ E. Using the focal points of E, which are located at one can introduce the two unit vectors v 1 and v 2 , We then consider the unit vector v which is a quadratic equation in r , The point p that we look for is the closest to c. Therefore, we choose the smallest of the two roots of the quadratic equation. Introducing the positive quantities α = (v T Dv), β = −(v T Dc), and γ = (c T Dc), we thus obtain so that Other weighted averages of v 1 and v 2 could be considered, but the one in (90) seems to work best for the majority of configurations. In fact, the normal along an ellipse at p bisects the angle f 1 p f 2 and the simple average (90) is meant to mimic this observation, although in our algorithm p is approximated by c. The optimal weighted average would actually be the one such that the normal vector to E at p share the same direction as v since v would then be the normal unit vector to the circle centered at c and of radius r , i.e., tangent to E. In that case, the point p would actually be the solution x to the minimum potential pair approach (see Lemma 4). Note that in the case of two arbitrary ellipses, one could use the mapping of Sect. 4.2, find the guess point p using the algorithm above, and map back p to the current configuration using transformation (88).
In the same configuration of E and C as before, Džiugys and Peters [12] also proposed conditions to find an initial guess point p; see

Contact detection algorithms
The main goal of this section is to review the main contact detection algorithms for pairs of ellipses that have been proposed in the literature and, more specifically, recast the methods as minimization problems such as those satisfied by the MDP and MPP introduced in Sect. 3. Many of the algorithms to be discussed were not explicitly defined as minimization problems (with or without explicit constraints) and were often categorized differently by the researchers Fig. 12 The initial guess point p = ( p x , p y ) suggested by Džiugys and Peters [12] can be chosen as any arbitrary point satisfying the conditions (92) Fig. 13 The points x i and x j are the intersection points of ellipses E i and E j , i.e., We observe that the contact point x c , obtained here by the Intersection Algorithm (IA), is in this case close to the co-gradient locus H i j themselves. Since it is ultimately the authors' hope to better highlight the similarities and differences between the published algorithms, it is incumbent on us to introduce a new classification which may conflict with those found in the literature. Whenever possible, we will indicate those conflicts in naming and justify the new terms. The framework will consist of a pair of ellipses E i , E j ⊂ R 2 in near-perfect contact, with or without overlap, as defined in Definition 6, but we shall discuss, when deemed necessary, the behavior of the algorithms for other configurations of the ellipses. A common feature of all algorithms presented here is that they compute two points x i and x j on the ellipses E i and E j , respectively. Then, the contact point is usually defined as the midpoint between x i and x j and allows one to compute a penetration (or separation) distance d i j = x i −x j . Moreover, one can identify the contact normal as the unit vector from the normal vector n i (x i ) and the opposite normal vector n j (x j ) Alternatively, the contact normal could be computed as Another definition for the normal vector [34] has been provided as the unit vector along the line passing through the centers of the two tangent circles at x i and x j to the ellipses E i and E j , respectively. In the case of the intersection of two ellipses, the normal at the point of contact is usually defined as the unit vector perpendicular to the line passing through the two points of the intersection set [42]. Finally, the tangent line is defined as the line perpendicular to n c . As the relative positions of the ellipses approach that of perfect contact, then the points x i and x j coalesce, the distance d i j vanishes, and the tangent line aligns with the tangent to E i at x i and the tangent to E j at x j .

Intersection algorithm (IA)
The intersection algorithm was introduced by Rothenburg et al. in [42] and relies on estimating the points in the intersection set I i j = E i ∩ E j introduced in Definition 2. It is perhaps the most intuitive algorithm. We thus believe that it needs to be included in this survey, despite some of its drawbacks and limitations described below.
Lemma 2 shows that the intersection algorithm can be cast as the minimization Problem (39), which in the presence of overlap is the same as the minimum distance pair Problem (42) but without a constraint on the normals: From a practical point of view, the method, in either case, consists in solving the system of equations: where f i and f j from Eq. (13) are the global geometric potentials of E i and E j , respectively. A priori, we make no assumptions concerning the two ellipses. If E i and E j do not coincide, the system of equations can naturally be reduced into a single quartic equation in the first coordinate x (alternatively, in the second coordinate y): where the coefficients a k , k = 0, . . . , 4, of the polynomial can be explicitly given in terms of the coefficients of (96), see [42]. The quartic Eq. (97) admits at most four real roots x , = 1, . . . , 4. Moreover, for each root x , one can use an explicit formula, see [42], to compute the corresponding coordinate y .
Different configurations of the ellipses will provide a different number of solutions to the quartic Eq. (97): 1. There are no real root. This is the case when the two ellipses are disjoint, i.e., I i j = ∅, whether they are disjoint with non-penetrating CoM or one ellipse lies inside the other (penetrating CoM). In this case, the algorithm does not provide information about the penetration distance. 2. There is one real root of multiplicity two. This is the case when there is a point at which the two ellipses are in perfect contact or when the x-coordinate of two distinct intersection points coincide. In the latter case, one cannot compute the y-coordinate of the two distinct points using the explicit formula [42]. Instead, given the common coordinate x , one should solve the quadratic equation in y using the global potential of ellipse E i , Alternatively, one could consider solving the quartic equation in y to obtain two distinct roots y . If a single root in x corresponds to a single point (x, y), then the separation distance and the normal can be computed easily. 3. There are two real distinct roots only. This is the most common configuration if the ellipses are in near-perfect contact with overlap. If x i and x j are two intersection points, the length x i − x j provides an approximation of the length of the overlap. However, more work is needed to provide an estimate of the penetration, see more details in [42]. This case is illustrated in Fig. 13. 4. The four roots of the polynomial are all real, either all distinct, or two distinct roots and one root of multiplicity two, or two roots of multiplicity two. In practice, these cases are unlikely to occur in DEM applications since all are indicative of relatively large overlaps between ellipses.
Drawbacks in the intersection algorithm are that one needs to compute all real roots of the quartic polynomial and handle all specific cases depending on the number of roots and their multiplicity. Moreover, a major issue is that the algorithm is prone to numerical instabilities [42] in the case of very small overlaps between ellipses when estimating the two points {x i , x j } and the normal direction. For these reasons, it is usually not the recommended approach for contact detection. Finally, the extension of the intersection algorithm to 3-D is increasingly more elaborate since the intersection set I i j consists of 2-D ellipses on the surface of the ellipsoids. Given the weaknesses of this approach, we will not be providing further details on its extension to 3-D. However, we refer the reader to Ouadfel and Rothenburg [38] who proposed an algorithm in 3-D to find the contact point and contact normal.

Geometric potential algorithms (GPA)
A geometric potential algorithm was first described by Ng et al. [31,35] and has been further improved by Mustoe and Miyata [34]. GPAs are based on the symmetric pair of minimization problems (48)-(49) that we recall here for convenience For any pair of ellipses in near-perfect contact the two problems have unique solutions, see Lemma 4. As the numerical experiments will show, GPAs are quite robust, even for pairs of ellipses with high aspect ratios, and at a computational cost competitive with other methods. However, two distinct problems must be solved and each problem generates up to four critical points from which the global minimum must be found. We present several numerical implementations of the GPA, including a Lagrangian and a parametric formulation. Penalization has not been found to be an effective means of solving (48)-(49) because the wide range of values of volume and aspect ratio make the choice of the penalization parameter difficult. In GPA algorithms, x i and x j are found by solving two distinct problems in the same way. Therefore, in the following, we discuss only the problem of finding the point x j . The point x i is found in a similar fashion (Fig. 14).

L-GPA: The Lagrangian approach
We describe the solution method proposed in [31] for solving the constrained minimization problem (99) to find x j . We thus introduce the Lagrangian where the constraint x ∈ E j , i.e., f j (x) = 0, is enforced via the Lagrange multiplier λ. Using the representation (1) for f k , k = i, j, i.e., the stationary points (x, λ) of L i satisfy the system: Isolating x as a function of λ from (100), we find Substituting (102) for x(λ) in (101) produces a quartic polynomial in λ: whose coefficients a k can be computed explicitly and can be found in [31]. Solving (103) provides at least two and at most four real roots λ , which in turn yields four candidate points x(λ ), = 1, . . . , 4, according to (102). Then, the point x j is selected as the point that minimizes f i (x) among the points x(λ ). Extension of the Lagrangian approach to the case of ellipsoids is straightforward and results in a root-finding problem equivalent to (103), but involving a polynomial of degree 6 (see [31] for details).

P-GPA: The parametric approach
We briefly describe here the approach proposed in [34] to solve the minimization problem (99). The main idea is to use the parametric representation (15) of an ellipse in order to eliminate the constraint from the minimization problem. This technique works in both two and three dimensions. Let the local potential f i of E i be given as in (5), i.e., where The points x on E j can be parameterized in the local coordinate system by ζ (t) = (a j cos t, b j sin t) for t ∈ [−π, π[. Using (4), the coordinates of x in the local reference system (O, ξ, η) associated with E i are then given by Replacing ξ (t) in (104) by (105), one can then express f i as a function of parameter t only, i.e., which can be reduced to Fig. 14 The points x i and x j obtained by geometric potential algorithms that provide the MPP, i.e., (x i , x j ). Note that the contact point x c does not necessary belong to the gradient locus H i j Fig. 15 Graphs of the functions f i (t) (left) and f j (t) (right) associated with the ellipses E i and E j shown in Fig. 16. The function f i (t) has exactly one global minimum and one global maximum while the func-tion f j (t) has two local minima and two local maxima in [−π, π[. In both graphs, the local minima and maxima are represented by a dot (•) and a disc (•), respectively where θ i j = θ j − θ i , with θ k for k = i, j, being the angle of rotation of R k , as defined in (2). It follows that the constrained minimization problem (99) for x j can be recast into the unconstrained minimization problem in one variable from which one would obtain the point x j ∈ E j by the change of variable (4): The nonlinear functions f k (t), k = i, j, may have several extrema, which makes root-finding algorithms for ( f k ) (t) = 0 difficult. For example, in Fig. 15, we plot f i (t) and f j (t) associated with the pair of ellipses E i and ellipse E j shown in Fig. 16. In particular, we observe that the function f j (t) has two local minima and two local maxima. In order words, without any additional con- Fig. 16 Configuration of the two ellipses E i and E j for which f j (t) has multiple local minima, as shown in Fig. 15. The points corresponding to local minima and local maxima are represented by a dot (•) and a disc (•), respectively straint, one needs to search for all extrema in order to find the global minimum. 3. M-GPA: The mapping approach Džiugys and Peters [12] proposed an alternative approach by introducing the mapping described in Sect. 4.2, which transforms E i into a unit circleĈ i and the other ellipse E j into an ellipseÊ j in its local reference system by the same mapping. Assuming the map has been applied and dropping the hat symbols, the circle and ellipse are then given by We first observe that the E i -norm is simply the Euclidean norm as the ellipse E i is now reduced to the circle C i ; see Lemma 13. It follows that the minimization problem (99) now reads In other words, the problem is to find the point x j on ellipse E j that is the closest to c i in Euclidean norm. Džiugys and Peters [12,13] proposed to solve the minimization problem (108) using two different approaches, one based on an iteration method and the other based on partial analytical results. We shall present only the latter approach below. Introducing the distance ρ i from the center of circle C i to any point they explicitly enforce the constraint x ∈ E j by rewriting ρ i as a function of x only, using the equation of ellipse E j where κ = b j /a j and c i = (c x , c y ). The goal is therefore to find the global minimizer of the functional ρ 2 i , i.e., Finding the critical points of the minimization problem leads to solving a quartic equation in x: The quartic equation has a maximum of four roots, with possibly some conjugate complex roots, which can be found using an iterative nonlinear solver. Džiugys and Peters [12] also suggested an approach in which the solution procedure could be reduced to solving a cubic equation after introducing a special change of variable.
Once the roots x , for = 1, . . . , 4, are found, one can compute the corresponding y coordinate, and the point x j is selected among the solutions x = (x , y ) such that it minimizes ρ i (x) in (109).

The constrained geometric potential algorithm (C-GPA)
We have seen in Lemma 4 that the additional constraints on the normals (50) and (51) are always satisfied at the minimum potential pair (MPP). Although the additional constraints on the normals are non-binding and could simply be ignored, it does significantly modify how the problems are solved and should further stabilize the problems. Ting and his collaborators [44][45][46] proposed to relax the constraints on the normals and to simply enforce that the points of the MPP belong to the so-called co-gradient locus. The constrained geometric potential algorithm can therefore be viewed as an extension of the GPAs. In other words, the point x j (resp. x i ) of the MPP is defined as the closest point with respect to the E inorm (resp. E j -norm) to the center of the ellipse E i (resp. E j ) that intersects the co-gradient locus H i j and the ellipse E j (resp. E i ). Formally, the problems can be recast as the constrained minimization problems: x j = argmin We emphasize here that above problems are different from the minimization Problems (53) and (54). Indeed, the fact that a point x belongs to the co-gradient locus H i j does not necessarily imply that n i (x) + n j (x) = 0 as one could also have n i (x) = n j (x). The method proposed by Ting et al. in [45,46] to solve problems (111) and (112) aims at finding the intersection points between each ellipse and the locus H i j (the authors refer to H i j as the locus of common slope in two dimensions). Since the two branches of the hyperbola may cross twice each ellipse, the solutions x i and x j are chosen as those that minimize f i and f j (equivalently, the E i -and E j -norms), respectively. In order to simplify the analysis to search for x i , Ting et al. in [45,46] proposed to consider the transformation of Sect. 4.1 that maps E i into the unit circle C i centered at the origin. Similarly, the second point x j is found by transforming the other ellipse E j into the unit circle C j while mapping the ellipse E i into a new ellipse under the same transformation, see Fig. 17.
We briefly describe the algorithm in the case of a unit circle C j and an ellipse E i given by The equation of the co-gradient locus H i j , see Eq. (68), is given in that case by Fig. 17 The two steps of the C-GPA algorithm are illustrated above. First, the mapping of Sect. 4.1 is applied to transform ellipses E i and E j to ellipse E i and the unit circle C j with center at the origin. Second, the solution of Problem (112) provides the point x j as the closest point to ellipse E i among the points that intersect co-gradient locus H i j and circle C j where A is the anti-symmetric matrix (67). Using the notation of Sect. 2, the equation can be explicitly written as The intersection of H i j and C j can be found by combining this equation with x 2 + y 2 − 1 = 0 to obtain a single quartic Solving for (114) yields a maximum of four real roots x , = 1, . . . , 4. Ting et al. [45,46] compute the second coordinate y from x using the formula derived from the equations of the circle and ellipse: as long as the denominator does not vanish for x . There may be an issue when the solution of (114) yields one real root of multiplicity two, which is the case when the two intersection points between the circle and the ellipse have the same coordinate x . Similarly to IA, the computation of the ycoordinate using (115) would result in one intersection point. One remedy is to use the equation of the circle to solve for y y = ± 1 − x 2 .

The steered geometric potential algorithm (S-GPA)
The authors proposed a novel approach [29] for contact detection between pairs of ellipses and ellipsoids combining existing techniques and new features. The algorithm belongs to the class of geometrical potential methods and, more specifically, involves solving the constrained minimization problem (112) after applying the mapping of Sect. 4.2.
The efficiency of the algorithm thus relies on the following key ingredients: 1) the transformation (88) that maps the pair of ellipses into an ellipse centered at the origin and a unit circle; 2) the construction of an effective initial guess following the approach described in Sect. 4.3; 3) the use of Newton's method for the root finding problem; and 4) the introduction of an additional constraint to guarantee convergence to the desired root. Suppose now that the pair of ellipses, after transformation, consists of the ellipse E j centered at the origin in its local coordinate system and of the unit circle C i centered at c i = (c x , c y ). We are thus looking for the point x j that satisfies: Moreover, x ∈ E j can be expressed in parametric form with t ∈ [−π, π[ as Therefore, points x belong to E j ∩ H i j if and only if they are the roots of the nonlinear function in t: We first note that the roots of h(t) are actually the critical points of the function f i (x) when x ∈ E j . Indeed, introducing the function g(t) defined for t ∈ [−π, π[ as it is straightforward to show that h(t) is associated with the derivative of g(t) as follows: In other words, the roots of h(t) also satisfy g (t) = 0, meaning that the points in E j ∩ H i j are either at a local minimum or at a local maximum of distance from c i when travelling along the ellipse E i . Since we have seen that the hyperbola H i j may intersect ellipse E j up to four times, we know that there exist up to four roots of h(t). The point that corresponds to the global minimum of g(t) , and hence f i (x(t)), is thus the point x j = x(t j ) that we are looking for. The scalar function h(t) is clearly continuous on t ∈ [−π, π[ but obviously non-convex. If one wants to use the second-order Newton's method to find the roots of h(t), one needs to define an accurate initial guess point and possibly an additional constraint to ensure that the method converges to the desired root t j , which will provide the actual solution x j to Problem (112). The initial guess point provided by the approach described in Sect. 4.3, although generally accurate, is not guaranteed to be located inside the basin of attraction of t j , in which case the method will converge to another root of h(t).
Among the roots of the function h(t), only t j satisfies n i (x(t j )) · n j (x(t j )) < 0 or, in an equivalent manner, ∇ f i (x(t j )) · ∇ f j (x(t j )) < 0; see [29] for a formal justification. This motivates us to introduce the function q(t), defined on [−π, π[, as and the constraint set Therefore, the solution x j to Problem (112) can be obtained by searching for the unique root t j of h(t) that satisfies the constraint t j ∈ S. Using Newton's method, we check that each new iterate belongs to S. If not, we apply one iteration of the line search method to get closer to t j , and we assume back in S, hence the name chosen for this approach. Extension of the algorithm to the case of pairs of ellipsoids is straightforward and is described in [29]. The advantage of this approach is that only the global solution to the minimization problem is actually computed and the use of the Newton's method and the choice of the initial guess point make it very efficient, as shown in the numerical examples of Sect. 6.5.

The common normal algorithm (CNA)
The common normal algorithm (CNA), first introduced by Lin et al. [31] in the case of ellipsoids, aims at rewriting the condition on the normals (43) satisfied by the MDP (42) into a system of solvable equations for the points x i ∈ E i and x j ∈ E j . We recall that the condition on the normals, namely is a natural property of the minimum distance pair of closest points (x i , x j ) ∈ E i × E j when the ellipses are disjoint, but the constraint is actually necessary when the two ellipses overlap, as the minimization problem (42) without the constraint on the normals leads in this case to the intersection set. Our objective in this section is to re-interpret the original formulation of the CNA [31] as a minimization problem. Let E i and E j be two arbitrary ellipses. For any given point x i on E i , one can always identify a unique point x j on E j such that n j (x j ) + n i (x i ) = 0. This implies that the set of constraints (x i , x j ) ∈ E i × E j and n i (x i ) + n j (x j ) = 0 provides an underdetermined system and, more precisely, one additional scalar constraint is required. Lin and Ng [31] proposed to consider that the unit vector going from x i toward x j be also equal to the normal vectors n j and −n i . Introducing the unit vector the problem of finding x i and x j would then consist in the following set of equations: Unfortunately, the problem above presents a few issues, which can be clearly described in the case of two circles: (1) If the two circles are disjoint, the problem has a unique solution pair, but the distance between the two points x i and x j reaches a maximum rather than a minimum, meaning that it is not really the pair of points that one is looking for; (2) if the two circles overlap, the problem admits two solutions, one solution pair for which the distance between the two points is a minimum and the other for which the distance is a maximum; (3) moreover, if the overlap between the two circles becomes very small, the calculation of the vector n(x i , x j ) defined in (118) becomes problematic and the problem is ill-posed in the limit case of perfect contact. We illustrate in Fig. 18 the existence of two solutions to the system of equations (119). Moreover, the system (119) consists of six nonlinear equations for the four variables x i = (x i , y i ) and x j = (x j , y j ). In order to circumvent the issue, Lin and Ng [31] proposed to consider only the following four equations (adapting the method described in 3-D in [31] to the 2-D case): Unfortunately, this arbitrary reduction in the number of equations has the effect of introducing additional solutions, as shown in Fig. 19. This is due to the fact that one should consider the direction of the normal vectors rather than simply equating their components in the x-direction, as done in (120). For all these reasons, the common normal algorithm is deemed inappropriate for finding contact points between ellipses or ellipsoids. As a final remark, the constrained minimization problem associated with Problem (119) can be cast as We note that the above problem is actually similar to the minimization problem (47) except for the choice of the objective function. We now elaborate on the closest co-normal algorithm.

The closest co-normal algorithm (CCA)
The algorithm was proposed in [53] as a variation of the common normal algorithm. The problem can be formulated as the closest co-normal minimization problem (47) that we recall here for convenience: We observe that the points x i and x j are not necessary located that is, so that Using the above formula, one can now express the difference In order to apply the constraint n i +n j = 0, one can introduce the common direction n, parameterized with respect to t ∈ [−π, π[, i.e., such that n(t) = n j = −n i . Therefore, The minimization problem (122) can thus be recast as that of finding the value t i j ∈ [−π, π[ such that The original minimization problem with multiple constraints in the four variables (x i , x j ) is now replaced by the unconstrained minimization problem (129) in the only variable t.
We show in Fig. 20 the solution pair (x i , x j ) obtained by solving (129) and observe that the two points do not necessarily belong to the co-gradient locus H i j . However, we note that for some configurations of ellipses, the distance d(t) = x j (t) − x i (t) may be difficult to minimize as it may exhibit very flat regions and multiple extrema (up to two local minima and two local maxima), as shown in Fig. 21.
Finally, we note that the algorithm extends straightforwardly to the case of ellipsoids by considering the following parameterization of the normal vector: (130)

Numerical results
The objective of this section is twofold. First, it is to provide examples where one can observe the significant differences  between the contact point associated with the intersection set (IS), the MDP, and the MPP. Second, it is to provide direct comparisons in accuracy and efficacy between the different contact detection algorithms presented in Sect. 5. In fact, the existing literature only provides a few comparisons between existing algorithms, mostly on a few pairs of ellipses, and those comparisons never emphasize how different the contact points from different algorithms may be. In practice, comparisons are difficult because (i) algorithms based on different minimization problems have fundamentally different solutions; (ii) some algorithms may be based on root finding of scalar polynomial equations, others on large coupled nonlinear systems, while others may utilize a key transformation to normalized coordinate systems. These choices dramatically affect accuracy and cost, how they are initialized, and even how accuracy is measured; (iii) the wide variety of numerical techniques deployed imply that the computational efficiency is highly dependent on specific implementation details (code optimization, language strengths, or platform choice); (iv) some algorithms have weaknesses in robustness, or may not be able to exploit prior accurate estimates of contacts, which render them undesirable in certain applications.
These difficulties will be circumvented by focusing mostly on the underlying minimization problems driving the algorithms and then dedicating the last subsection to applying the algorithms on a randomly generated set of pairs of ellipses in close or almost perfect contact. The first four examples, detailed in Sects. 6.1-6.4, will present pairs of ellipses for which the differences are attributable solely to the minimization problem on which they are based. This implies that the computational cost will be ignored and all the algorithms will be run with a high tolerance in order to provide to machine precision the exact solution to the minimization problem. For these first four test problems, none of the contact detection algorithms failed, so we could focus only on the minimization problem.
In Sect. 6.5, our objective is to highlight the algorithmic aspects of the resolution of the different minimization problems. To accomplish this, ten thousand pairs of ellipses in close or almost perfect contact are randomly generated according to an algorithm detailed in another publication by the authors [29]. For each algorithm, the total computational time is found and then further subdivided into its different steps; see Tables 9 and 10. This allows us to address issues (ii) and (iii). For those tests, two different tolerances were used. Overall, Sect. 6.5 indicates that the S-GPA is the fastest, but because of the issues just mentioned above any comparison between algorithms comes with significant caveats.

Comparison of the intersection set (IS), the MDP, and the MPP
The objective of this example is to demonstrate that with strict tolerances, the contact point obtained by a given algorithm will depend only on the minimization problem which the algorithm attempted to solve. Following this remark, the numerical results of Sects. 6.2-6.4 will focus only on the different solutions associated with either the MDP or the MPP. Nevertheless, the example in this section is also interesting in its own right since it presents a pair of overlapping ellipses, given in Table 1, for which the intersection set (x i , x j ), the minimum distance pair (z i , z j ), and the minimum potential pair ( y i , y j ) are different. In particular, this example shows that the resulting contact points x c , z c , and y c will be dis-tinct, as seen in Fig. 23. On the other hand, the normals at each of these three contact points are roughly the same, computed according to either (94) or (95). In general, we have observed the normals to be less sensitive to the choice of the algorithm than the estimates of the contact point. The numerical results in Table 2 show that the MPP ( y i , y j ) is the same whether it is computed by the P-GPA, L-GPA, M-GPA, or the C-GPA. Similarly, using sufficiently high tolerances, we find that the estimated MDP is the same whether computed by the CNA (when CNA provides the correct answer) or the CCA.
Examining Fig. 23 and comparing it to the definition of the MDP and MPP, it is easy to explain the differences between the contact points x c , z c , and y c . For example, the MPP ( y i , y j ) clearly belongs to the co-gradient locus, but their normals are different at each point. On the other hand, it is clear in Fig. 23 that the normals at the MDP (z i , z j ) are opposite. The contact point z c is far from the co-gradient locus because both ellipses have roughly parallel surfaces along a large portion of the intersection E i ∩ E j .

Different MDP and MPP with the same contact points
In this example, we present a pair of ellipses, described in Table 3, for which the contact point z c of the MDP (z i , z j ) and the contact point y c of the MPP ( y i , y j ) are the same, but the associated normals are different. This expands on the previous example because the pairs are different, yet lead to equal contacts. The fact that the normals are different but the contact points are the same implies that the choice of using the MDP or the MPP could lead to significantly different forces between ellipses in a DEM model. The example is quite simple because, as Table 3 shows, both ellipses have the same aspect ratio and the co-gradient locus degenerates to a straight line through both centers. The quartic equations (110) and (114) derived from M-GPA and C-GPA, respectively, degenerate in this case to quadratic equations, i.e., a 4 = a 3 = a 1 = 0. This implies that the mappings in Sects. 4.1 and 4.2, applied, respectively, in M-GPA and C-GPA, send both ellipses to circles. In those new coordinates, the contact is easy to compute analytically. Those analytic verifications show that the coincidence z c = y c is exact, and not simply an artifact of the numerical algorithms. We report in Table 4 and show in Fig. 24 the MDP and the MPP, obtained, respectively, by a GPA algorithm and by the CCA with high tolerance. It is a coincidence in this example that the normals are opposite, both for the MDP and for the MPP. Table 2 The contact points x c and their normal vectors n c for the different contact detection algorithms in Example 6.1 Algorithm x i x j x c n i n j n c Fig. 23 Location of the contact points using different contact detection algorithms: the intersection set is the pair (x i , x j ) with contact point x c , the MPP is ( y i , y j ) with contact point y c , and the MDP is The co-gradient locus H i j is drawn as a dashed line, and it traverses both centers c i and c j . The dotted line presents the scaled ellipses that are tangent at the MPP Table 3 The two ellipses E i and E j of Example 6.2

Ellipses in perfect contact
Following the work of Dziugys et al. [12,13], we consider a pair of ellipses of high aspect ratio in perfect contact, according to Definition 4. This is a numerically challenging case, yet both the MDP and the MPP coincide mathematically in this case. As expected, all algorithms identify the same contact point. This perfect contact is illustrated in Fig. 25. The ellipses are described in Table 5, and the resulting MDP and MPP pairs are given in Table 6.  , y j ). However, the contact points z c and y c are identical and are both located on the co-gradient locus H i j Fig. 25 The location of the contact points x c is the same for both the MDP and MPP pairs z c associated with the MDP belongs to the intersection of both ellipses. However, the MPP produces a reasonable contact point y c inside the intersection between both ellipses E i ∩ E j . We insist here that these ellipses are not in nearperfect contact, that is according to Definition 6, because the two disks of radius ρ i and ρ j , tangent to y i and y j , respectively, do not have non-penetrating CoM, as required by Theorem 7. The description of the ellipses is given in Table 7. The estimates of the contact points for both the IS, the MDP, and the MPP are given in Table 8. We note that the estimated normals are very close to each other.

Large-scale tests
In contrast to the previous four sections, we attempt here to analyze and compare the accuracy and efficiency of the algorithms presented in Sect. 5. More precisely, we study the individual numerical approximations implemented in the algorithms used to solve the two main contact detection problems we identified, that is the MDP and the MPP, independent of the characteristics of the problems themselves studied in the previous four sections. As we argued in the introductory comments of this section, there are many reasons why comparisons between contact detection algorithms could be difficult. Nevertheless, as we will explain below, we can come to a limited number of conclusions by studying the behavior of these algorithms on a large sample of pairs of ellipses in close contact, and then separately considering the contributions to accuracy and efficiency of the individual components of the algorithm. The tests in this section are new to the literature and should help to establish benchmarks for comparisons between such algorithms for contact detection. The first step was to generate a random set of 10,000 pairs of ellipses in almost or close contact, according to an algorithm described in a companion publication [29], and to apply to each pair of ellipses all algorithms (except IA) of Sect. 5. In order to reproduce pairs of ellipses one might encounter in DEM, the generating algorithm provided some control on the aspect ratio of each ellipse, on their relative orientation, on their relative closeness, and on the location of the contact point along the boundary of each ellipse. First of all, the first ellipse E i was permitted to have maximum aspect ratio of a/b = 5, while ellipse E j was permitted to have an aspect ratio as large as a/b = 20. Figure 27 presents the actual distribution of the aspect ratio a/b for E i and E j . The generating algorithm assumes that a j b j = 1 for E j , but the distribution of the area πab for E i is randomly determined and shown in Fig. 28.
The algorithm generating these pairs is able to provide exactly one of the two points in the MPP, thus allowing us to estimate the relative error of the algorithms in the GPA family. On the other hand, for the IA and the MDP algorithms, estimates of error are not directly accessible. In any case, we are able to control the distribution of the penetration/separation distance, as shown in Fig. 28, supporting our claim that the pairs of ellipses in our tests were relatively close. Furthermore, the generation of the MPP for the pair of ellipses also allows one to uniformly select the second ellipse along the boundary of the first ellipse, thereby ensuring that we consider contact points both in the flatter or in the more curved regions of the boundary of the ellipses. The numerical results of the experiments are summarized in Table 9, where a relative tolerance of 10 −5 was used, and in Table 10, where a stricter relative tolerance of 10 −9 was used. We will discuss the results of Table 10 later, since it mostly concerns the observed convergence and how it depends on the underlying numerical approximations. Table 9 presents for each of the seven algorithms of Sect. 5, the total computational time required for the resolution of the 10,000 pairs of ellipses, the statistics of the number of iterations required for the solution, and the statistics of the error in those approximate solutions. First of all, the error could only be measured for the algorithms estimating the MPP, since the algorithm generating the pairs only provided the exact MPP. This explains why the errors for the CNA and the CCA are not tabulated. The table also includes the number of iterations  ( y i , y j , y c ) are obtained by using GPA, and (z i , z j , z c ) are obtained by using CCA Table 8 Contact point x c and normal vector n c obtained by different algorithms in Example 6.4 Table 9 Computational cost, number of iterations, and logarithmic error for different algorithms using a relative tolerance of 10 that were required to attain the desired relative tolerance, but one must be careful when comparing these values because the nature of the iterations in, say, the golden search algorithm, Newton's method, or Francis' algorithm is completely different. Finally, we have included the total computational time required to solve the 10,000 pairs of ellipses, as measured by the profiler in MATLAB. Although estimates of computational time in MATLAB are known to be somewhat variable, we have performed many such studies and found the estimates of computational time to be consistently reproducible to within 10%. We now proceed to analyze the results of the experiments in Table 9, going from the most general to the most specific conclusions. First of all, we remark that the median error was roughly 10 −11 , that is several orders of magnitude lower than the chosen tolerance, showing that for most pairs of ellipses, the algorithms converged quickly to the MPP. The relatively low standard deviation further shows that the error was close to the relative tolerance of 10 −5 for only a small subset of the pairs of ellipses. In Table 9, we immediately remark that the total computational time is only roughly equal to the sum of the time required for the different components because we omitted the computationally insignificant but necessary step of computing the coefficients in the systems of equations we needed to solve. We observe that in the L-GPA, M-GPA and C-GPA algorithms, a significant fraction of the computational effort is spent on finding the roots of a polynomial. However, in the S-GPA algorithm, mapping costs more than finding a root by Newton's method. The data also show that the most expensive algorithm was CNA because it required a system of equations to be solved rather than simply finding the roots of a polynomial of degree four as in the case of the algorithms L-GPA, M-GPA, and C-GPA. Table 9 contains some interesting observations about the efficacy of the iterative solvers underlying some of these algorithms. First of all, it appears that the Lagrange approach leads to the system requiring the smallest number of iterations. On the other hand, the P-GPA uses a golden search algorithm to compute initial estimates of the minima and maxima before invoking Newton's method to converge rapidly to the minima and maxima. Unfortunately, our tests indicate that it is difficult to reduce the number of iterations in the golden search without obtaining initial estimates for which Newton's method will not converge. As a matter of fact, even using the golden search algorithm with the recommended tolerances, roughly 1% of the pairs of ellipses did not converge to a contact point for the P-GPA. In order to maintain the consistency of our tests, the pairs of ellipses for which P-GPA failed to converge were also excluded from our tests with the other methods. We remark that the standard deviation of number of iterations for M-GPA and C-GPA are higher than the other algorithms, which leads us to conclude that the number of iterations of M-GPA and C-GPA depend on the relative configuration of the two ellipses. In contrast, the convergence of the S-GPA appears to be independent of the relative geometry between the ellipses.
Later in this section, we will examine the pairs of ellipses for which different algorithms attained the maximum number of iterations; see Figs. 29, 30, and 31. Although the number of iterations required by different algorithms is not necessarily correlated, the pairs of ellipses for which a specific algorithm Overall, the data indicate that the S-GPA algorithm was the most efficient. Our hypothesis is that this algorithm combines a good initial estimate of the MPP (for an ellipse at origin and a unit circle) with Newton's method quadratic convergence. In practice, we found that the S-GPA often converged in a single iteration. Finally, the numerical experiments indicate that the CNA and CCA algorithms were by far the most costly alternatives. We will therefore refrain from discussing them any further.
We now turn to the data in Table 10 obtained using the same set of pairs of ellipses but with a stricter tolerance of 10 −9 . A priori, we expect the results to indicate the same overall trends but the stricter tolerance should help to identify any robustness issues. It is clear that the stricter tolerance produces more accurate estimates and requires a larger number of iterations. Yet, it is noticeable that the median number of iterations is almost identical, while the median error is roughly 10 −4 smaller. This indicates that for the majority of the test cases we considered, the algorithm was already within the asymptotic regime of convergence and in many cases one or no iteration was required to satisfy the desired tolerances.
Among the L-GPA, M-GPA, and C-GPA, the L-GPA algorithm is again the most efficient and accurate, but it is still an order of magnitude slower than the S-GPA. The low standard deviation of the number of iterations suggests that the S-GPA was also the most robust algorithm. It appears that the M-GPA and C-GPA algorithms both required more iterations of their root-finding algorithm, Francis' method, in order to obtain the MPP, particularly when contrasted with L-GPA.
We hypothesize that the mapping step, present in M-GPA and C-GPA but not in L-GPA, might make the root-finding problem harder, although further tests would be required to confirm this.
Finally, we conclude this section by examining pairs of ellipses for which certain algorithms required the largest number of iterations from among our sample set. Figure 29 presents the pair of ellipses that required the largest number of iterations in Francis' method. In this case, the two ellipses appear to have their principal axis roughly aligned and to be in contact near the regions of the highest curvature. We also note that in this configuration, the points on the segment formed by y i and y j cross the segment formed of z i and z j . It is somewhat surprising that this configuration might be difficult for the M-GPA to handle. Figure 30 presents the pair of ellipses that required the most iterations of the L-GPA. This configuration appears very similar to the one associated with the M-GPA. Finally, Fig. 31 presents the worst-case scenario for the C-GPA and again the configuration seems unexceptional although the principal axis is angled by roughly π/4 and the contact occurs at points where the curvature is intermediate. In other words, the pair of ellipses in Fig. 31 is completely the opposite to what we found in the previous two configurations. We did not present the pair of ellipses associated with the worst-case scenario for the S-GPA algorithm because there was very little variation on the number of iterations. Furthermore, the S-GPA algorithm profited from a unique initialization; hence, a comparison with the other GPA-type algorithms would be biased.

Conclusions and perspectives
This research attempted to review and restructure the problem of detecting and estimating contact detection between pairs of ellipses. The topic has been studied by numerous engineers with most of the algorithms resting on heuristics, but as a result, much of the mathematical theory had yet to be fleshed out. This research has filled some of the mathematical and numerical gaps, particularly in 2D, and this has led the authors to reclassify published algorithms in order to highlight their differences and similarities. Furthermore, the authors have proposed a methodology to compare fairly the different contact detection algorithms (CDAs), something that was lacking in the literature. Finally, the authors demonstrated the usefulness of their analysis by introducing a novel CDA combining old and new techniques.
Although this research has identified and characterized several concepts that were implicit or ignored in past In this configuration of ellipses, C-GPA requires its largest number of iterations to find y i research, it is the authors hope that this work will help to stimulate research on this topic. Much of the mathematical theory still needs to be extended to 3D, particularly the co-gradient locus, and a proper study may provide some new insight into computational aspects of the subject. A careful analysis of the mathematical theory indicates that the main tools are the convexity of the ellipses and scaling arguments in projective space, which suggests that the notion of the co-gradient locus, and many of its associated concepts, could be extended to other families of convex particles. There are numerous mechanical, and even esoteric, procedures to construct ellipses, some of which may lead to original CDAs. This research only briefly mentioned connections to the Perram-Wertheim contact theory, or to the time-dependent contact detection problem, and connections to both of these theories could be fruitful. On the other hand, there exists a large literature on contact potentials [16,39] that may be very instructive to designing better CDAs, if not just better initial approximations for contact points. Our research has identified three notions of contact that could be used in collision detection, and although we believe these may be the only ones possible, there may exist others. This paper has discussed only the narrow phase of contact detection, but reductions in the computational cost of the broad-phase search could have significant impacts on DEM/MD simulations. It seems reasonable to assume that machine learning techniques could lead to breakthrough on this problem.
The numerical methodology for the comparison of CDAs could be pursued even further to obtain more precise information about accuracy and robustness. Although the notions of contact all coincide, it is not clear whether or not these notions might introduce long-term biases in large-scale particle simulations. On the other hand, the potential for bias is more likely in approximations of particles by clusters of spheres, and even that issue has not been explored. Finally, there is still a need for more test problems, such as those of Sects. 6.1-6.4, and we encourage other researchers to continue to share challenging configurations of ellipses. tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.