The Twinning Problem

In the field of crystallography, some crystals are not made of a single component but are instead twinned.In these cases, the observed intensities at some points in the lattice will be far larger than predictions. If we find the rotation associated to the twinned component, we can model this twin and improve our agreement with observations. In this report, we explore many routes to improve the process of identifying twins: Generation of fake data for better understanding and accurate testing. The representation of a rotation as defined by an axis and angle. The representation of a rotation as a quaternion. Using lattice points which must be equidistant from the origin to create our viable rotations. An algorithm focused on restricted possibilities. An exploration of 2D lattices for which twinning is mathematically impossible. We find that there is much to be investigated in the field of twinning.


Summary
In the field of crystallography, some crystals are not made of a single component but are instead twinned.In these cases, the observed intensities at some points in the lattice will be far larger than predictions.If we find the rotation associated to the twinned component, we can model this twin and improve our agreement with observations.In this report, we explore many routes to improve the process of identifying twins: • Generation of fake data for better understanding and accurate testing.
• The representation of a rotation as defined by an axis and angle.
• The representation of a rotation as a quaternion.
• Using lattice points which must be equidistant from the origin to create our viable rotations.• An algorithm focused on restricted possibilities.
• An exploration of 2D lattices for which twinning is mathematically impossible.
We find that there is much to be investigated in the field of twinning.

Introduction
OlexSys creates and maintains crystallographic software, enabling the structure solution and refinement of crystalline structures, through the use of X-ray diffraction.The problem presented by OlexSys to the study group concerns the post-hoc detection of twinning in the crystal diffraction pattern.Let us briefly explain the context of this problem.A crystal diffracts an incoming X-ray beam which leads to constructive interference of the outgoing diffracted beam.This interference yields particular peaks of intensities (so called Bragg peaks) which are arranged in a lattice structure.These peaks are detected by the diffractometer (Figure 1) and are available in form of given intensities at each of these lattice points.In other words, we have a lattice Γ * = Z a * + Z b * + Z c * with associated intensities I(h, k, l) at the lattice point h a * + k b * + l c * .In fact, these intensities are related to the internal structure of the crystal (its electron density) via a Fourier transformation.The challenge of structure determination is to derive from this information the internal structure of the crystal -that is, the kind of atoms and their positions in the unit cell of the crystal.This task is achieved by first obtaining an initial model, and then refining it with the software Olex2-refine.
In the real world, crystals sometimes are overlaid by two identical components which are related by a rotation.This phenomenon is called twinning.Due to the fact that both components contribute to the diffraction pattern, some of the intensities at certain lattice points are stronger than they would have been in the case of a single component.This happens precisely at those Bragg Peaks of both components which coincide.For that reason, the problem of detecting a twin can be reduced to a purely mathematical problem about rotations of the lattice Γ * .Namely, given a lattice Γ * , and potentially a subset of lattice points H with underestimated theoretical intensities, we wish to find rotation matrices R for which some points of RΓ * overlap, or are close to, points of Γ * , especially the points in H.These rotation matrices are called twin laws.An example of a twinned crystal can be seen in Figure 2, illustrating how the individual components may be interlaced.In the study group, we first focused on obtaining a deeper understanding of the problem at hand, with the creation of fake data (see Section 2) and observing the impact of twins on that data.It turned out in due course that rotations in 2D are much easier to handle than 3D rotations, as 2D rotations are uniquely determined by the rotation angle and the center.
In 3D, there is a 1-dimensional set of rotations mapping a point into another point at the same distance to the center.This fact makes the treatment of the 3D case more difficult.In our geometric investigations, we developed connections to quaternions and an interesting description of these sets of rotations by ellipses (see Section 3).
We investigated whether there were lattices Γ * in 2 dimensions for which non-trivial rotations have no overlaps, except for the origin.By non-trivial, we mean that we exclude the 180 • rotation, where every point is mapped to its negative.These investigations are presented in Section 5.As an intriguing byproduct, we found a connection between this and geodesics in hyperbolic geometry.
In this study group, we did not complete the development of an algorithm for finding twin laws in the 3D case, but the ideas which emerged provide various possibilities for such an algorithm.The challenge for the future is to develop such an algorithm which is particularly efficient.

Generating realistic fake data
In order to develop algorithms to detect twinning, an algorithm to develop realistic fake experimental data was developed.This has the advantage of being able to be quickly created and tweaked to the purposes of the study, but does not take into account any of the chemistry that would determine the results of real experiments.Both two-dimensional and three-dimensional algorithms were developed, and we present below in detail the algorithm for the two-dimensional case and subsequently some example results for both the two-dimensional and three-dimensional cases.

Two-dimensional fake data
We start by defining a 2D lattice with basis vectors a * , b * .We then generate a set of integer points {(h, k) ∈ Z 2 }, such that h a * + k b * represents a location in reciprocal space.For each (h, k), we generate a randomised weights w h,k = σ h,k exp(−(h 2 + k 2 )/α), where σ h,k is a uniformly distributed random number between 0 and 1, and α is a scaling factor.These weights determine the chance a X-ray event is detected at a point, and they decay for larger values of h and k, as seen in experimental data.
The fake data is generated by selecting a large number of points, N , to act as X-ray detection events.For each event, an (h, k) pair is selected uniformly at random using the weights w h,k .If we define the centre by (x * , y * ) = h a * + k b * , then the point is recorded at (x * r , y * r ) = N ((x * , y * ), Σ), where N is the multivariate normal distribution with covariance matrix Σ.The covariance matrix Σ is chosen to be a multiple of the identity matrix with small diagonal values so that most points fall within a small radius of the centre (x * , y * ).Once this procedure has been repeated for all N points, we are left with a set of randomised points that are clustered near the lattice points defined by h a * + k b * .
The intensity measured at each point (h, k) is determined by summing up the number of points that are measured within a radius R tol of the centre points (x * , y * ) = h a * + k b * .The result is a list of (h, k) points with associated intensities I h,k .
We can readily adjust this algorithm to observe twinning behaviour.Say we have a twinned component of angle θ with twin intensity β.We consider two sets of weights; the first, {w h,k , (h, k) ∈ Z 2 }, is generated as before and represent the main lattice components.The second, { wh,k = βw h,k , (h, k) ∈ Z 2 }, represent the twinned components.Now, when generating the N fake data points, the (h, k) pair is now selected using both sets of weights.This means for small values of β, most points will still be selected from the main component, but a certain selected percentage will be twins.If a twinned point is selected, the centre is given by rotating the point h a * + k b * by the twinning angle, and the point is generated close to this rotated location.
Example outputs of this algorithm are shown in Figure 3.The output is represented as a lattice of points, with their sizes proportional to the measured intensity.Figure 3a shows an example output for a square lattice with no twinning present.We can introduce a twinned component, at an angle 36.9°withβ = 0.1, and show the resulting intensities in Figure 3b.It can be seen that certain points now have higher values of the intensity, as a result of the overlapping behaviour of twinning.To show this, in Figure 3c we show in blue the original lattice rotated by the twinning angle.The green points on this figure are the difference in the intensities between the twinned and no-twinning cases.As can be seen, these points overlap with the rotated lattice, showing the effects of the twinning.
These examples show how it is possible to find the twinning angle by rotating the original lattice, and identifying the angles where many of the points overlap, as shown in Figure 3c.This can be posed as an optimisation problem, as we want to find the angle which maximises the number of overlapping points.Unfortunately, this is a difficult to optimise problem, as if we rotated the lattice slightly away from the twinning angle, the number of overlapping points drops to zero, so the search space of angles would have to be very fine.This method would also translate poorly to the three-dimensional problem, as the extra dimension in the search space would make this optimisation approach computationally prohibitive.

Three-dimensional fake data
In practice, the experimental data is three-dimensional.The algorithm for generating three-dimensional data is exactly the same as the two-dimensional case, except now considering (h, k, l) triples instead of (h, k) pairs, so for brevity the details are not included.As mentioned above, it is prohibitively expensive to brute-force numerically find the rotation matrix in three-dimensions, however this fake data is useful for visualising the problem and to make test cases for use with the mathematical models.Figure 4 shows an example of three-dimensional fake data on a cubic lattice.This time, all points are shown as the same size, regardless of intensity, for visual clarity.In three-dimensions, it can be difficult to see the lattice structure and the data can look disordered, as in Figure 4a.However the crystal structure can emerge by rotating the azimuthal angle, as seen in Figure 4b.The lattice rotated by the twinning angle, with the green dots showing the residual intensities between the twinned and no-twinning casces.
3 Rotations between two unit vectors in the 2-sphere Through this section, we will be describing rotations in 3D, how to represent them as quaternions, and how one can obtain a rotation from the knowledge that two points are mapped to one another.Any rotation centered on the origin (which is always the case for twin laws) can be written as a matrix R ∈ SO(n) such that R u 1 = u 2 (where R rotates u 1 onto u 2 ).It is also the case that R n = n, for n the 'axis of rotation'.

Rotations in 3D
We consider a lattice Γ = α a + β b + γ c, α, β, γ ∈ Z which is generated by three linearly independent basis vectors a, b, c.Two lattice points P 1 , P 2 ∈ Γ can only be rotated into each other if and only if both points have the same distance from the origin.It is assumed that the axis of rotation must go through the origin.Let p 1 and p 2 be the position vectors of the two lattice points, this means | p 1 | = | p 2 |.We also introduce the standard notation v for a unit vector.
For the following derivation, we require a definition of the cross-product matrix.Given a vector k, its cross product matrix is Given a unit axis of rotation n and a rotation angle θ ∈ [0, 2π) by the right hand rule, we can obtain the rotation matrix via Rodrigues' Rotation Formula: We now consider rotations which map p1 into p2 , note that p1 , p2 ∈ S 2 .The geodesic rotation gives the smallest possible angle between the two points, α = arccos(p 1 • p2 ), and has an axis of n1 = (p 1 × p2 )/ sin α.On the other end of the spectrum is the half-turn, a rotation with an angle of π and an axis of n2 = (p 1 + p2 )/( 2(1 + cos α)).This means any two points can be connected by many different rotations, if we do not insist on the Laura Midgley et.al.
geodesic rotation.The rotation angle will be bounded from below by the geodesic angle and from above by a half-turn.
Any permissible rotation axis is thus a linear combination of the axes n0 and n1 which we write now write as This vector is normalised by construction and depends on one parameter t ∈ [0, π).We describe the corresponding angle as θ(t), and now derive its dependence on t.
We begin with the knowledge that p1 • p2 = cos α, and use that p2 = Rp 1 .Thus Then, we note that the second term is zero, as (n(t)× p1 ) provides a vector perpendicular to p1 , and thus the dot product with p1 is zero.As p1 is of unit length, the first term is 1.Finally, we also make use of the fact that, given a particular scalar triple product a to reduce the third term which leads to We can further simplify the right hand side by substituting in n(t) = cos(t)n 1 + sin(t)n 2 .
Then we can distribute the cross product to obtain and as n2 is parallel to (p 1 + p2 ), p1 × n2 is parallel to n1 , and thus the two terms within this norm are perpendicular.This allows us to utilise Pythagoras' Theorem, and separate this norm: Plugging this into (3.4)yields
The verification of this relies on the specific tangent half-angle formula tan(θ/2) = 1−cos θ 1+cos θ .Then one can simply expand to verify: The positive and negative components represent the fact that we can take positive or negative theta to represent a rotation from p 1 to p 2 , or from p 2 to p 1 .In the initial equation (3.6), this could be chosen via taking positive or negative arccos.

Quaternion representations
Quaternions can been seen as extension of the complex numbers, we write a quaternion as q = q 0 + q 1 i + q 2 j + q 3 k where the unit vectors i, j, k satisfy The q 0 part of the quaternion is generally referred to as the real part, similar to complex numbers.It is well known that unit quaternions are useful quantities when representing spatial rotations in 3D.One, somewhat hand waving, way is to see this is as follows: Let (w, x, y, z) ∈ R 4 be a point in four-dimensional Euclidean space.Assuming the condition w 2 + x 2 + y 2 + z 2 = 1 makes this point an element of S 3 .A unit quaternion, an element of S 3 and an orthogonal matrix in 3D all have three degrees of freedom.Recall that an orthogonal matrix can be written as R = exp A where A is a skew-symmetric matrix, which in 3D has 3 independent components.
If n is a unit vector defining the axis of rotation and θ is the rotation angle, then the unit quaternion associated with this rotation is given by This formula is closely related to Euler's formula for complex numbers.Using the normal vector (3.3) and the angle (3.6) one would arrive at the canonical representation of our rotations.
If we consider an arbitrary vector v then its usual rotation would be given by v → R v for an orthogonal matrix R. We let q be the corresponding unit quaternion associated with R. Let us now map the components of v onto a quaternion w with vanishing real part q 0 = 0, then we would have w = v 1 i + v 2 j + v 3 k.The rotation given by the unit quaternion (3.8) now rotates our vector as follows: w → qwq −1 .This means rotations are realised as conjugations of w with q.
If we represent the imaginary part of our rotation q as a non-unit vector u, we can directly deduce that n = u/| u|, that q 0 = 1 − | u| 2 and that sin θ 2 = | u|.Such a non-unit vector u fully describes our rotation.A rotation with unit axis n and rotation angle θ is represented by u = sin θ 2 n.We now claim that the vectors u(t) = sin θ(t) 2 n(t) (with θ(t) and n(t) from the previous section) lie on the ellipse (3.9) We begin by noting that n(t) = cos(t)n 1 + sin(t)n 2 as given in (3.3), and hence u(t) is of the form u(t) = an 1 + bn 2 where we choose Taking a and b as given above and using sin 2 θ(t) 2 = 1 2 (1 − cos θ(t)), we obtain Thus we have determined that the set of rotations mapping p1 to p2 is represented by an ellipse with semi-minor axis sin α 2 n1 and semi-major axis n2 .We believe that this formulation could lead to further deductions on methods to determine twin laws.

Methods for Twin Finding
The study group has explored a systematic approach to try and identify the twinned structure of the reciprocal lattice by uncovering all possible rotation of this lattice that result in an overlapping of one or more lattice points.
In the literature, there exists an algorithm which iterates over axes and angles to brute-force locate viable rotations, see [2].The study group was interested in alternative methods to locate these rotations.
The rational behind our investigation was that since we always find the reciprocal lattice together with the intensified lattice points, which are the result of the twinned structure, we are able to identify this structure by first finding all possible rotations that may cause lattice points overlap, and comparing these rotation with the phenomenological data we posses.
Finding all aforementioned rotation will be done by an algorithm, which has the following observations underlying it: • In order for two lattice points to overlap completely by a rotation R their distance from the origin must be the same.If we'd like to allow for "radial tolerance" we can allow the points to be in a certain (small) spherical ring.• Once we restrict ourselves to an appropriate sphere/spherical ring, any two points can be made to overlap by an appropriate rotation matrix.In 2−dimension there is only one such matrix, yet in 3−dimension one find a one parameter family of such matrices.• Using these matrices, one can in fact count the number of points that will overlapgiving us a better possibility to discover the actual twinned structure in the original crystal.
• One can also allow for "angular tolerance", meaning that the rotations that are found are allowed to be perturbed slightly, to achieve more realistic results4

The overall algorithm
We are now ready to present our algorithm.
Step 1 (Listing the radii): Create a list G d that assigns to any lattice point its distance to the origin.We will commonly call these distances "Radii".
Step 2 (Identification of potential overlaps with radial tolerance): For any possible radius from the list G d and any ≥ 0, which represents the radial tolerance of our investigation, create the list G r, of all elements of G d whose radius is between r − and r +5 .If #G r, > 1 then there are two lattice points which could overlap under a particular rotation.
Finding these rotations and counting the number of overlaps calls for different tools in 2 and 3 dimensions.In both cases, however we define the list G r, ,∠ to be the pairs of lattice points from G r, together with cosine of the angle between them.Step 3 in 2D (Identification of possible rotations): Denoting by R θ the unique clockwise rotation by an angle θ we see that any pair of points from G r, ,∠ with the same angle θ between them (measured via the cosine) will create a unique overlap of points, up to the appropriate radial tolerance.Thus, the number of overlaps we will find under R θ equals the number of pairs of points from G r, ,∠ with the same angle between them.
Step 3 in 3D (Identification of possible rotations): As in the 2−dimensional case, any pair of points from G r, ,∠ yields a possible rotation that creates an overlap of points.Unlike the 2−dimensional case, however, there isn't a unique rotation mapping two lattice points to one another, up to the appropriate radial tolerance, but a family of such rotations (see Figure 5), defined by axes on the great circle lying between the two points.This family is identified uniquely by points on an ellipse (3.9) which is intimately connected to the pair of lattice points.
At this point, we consider two options: • Simultaneous considerations of two pairs of points yields a potential axis of rotation for twin laws -however, the angles must agree.The computation of the common axis is described in Section 4.2 • Searching for points in 3D space which have a high number of ellipse points close enough (up to some tolerance).A point that lies close to n ellipses corresponds to a rotation that acts on n pairs of lattice points, as as such creates n overlap patterns.
Step 4 (Adding angular tolerance): Create a list O r, that includes all possible rotations found in Step 3 together with all the pairs of lattice points which overlaps due to these rotations.For any given δ ≥ 0, which represents the angular tolerance of our investigation, create the list O r, ,δ which bunch elements of O r, whose rotations differ by an angle of δ (be it polar, azimuthal or a combination of the two in the 3−dimensional case).
We would like to conclude with a couple of remarks: • Note that the list O r, ,δ , attained in Step 4, not only gives us the list of all possible rotations with an radial tolerance and δ angular tolerance under which we will find an overlap of lattice points, but also all the points which will overlap.It is this information that we compare with our phenomenological reading in order to try and identify possible twinned structures.An example of this output for a 2D lattice is shown in figure 6. • Many crystals have inherent symmetry that will result in an "obvious" rotation that results in complete overlap of the lattice points.These rotations need to be discarded from the set O r, ,δ if want to find a "true" twinned structure.

Step 3 in 3D: Rotations in 3D for two pairs of points
We now look specifically at the finding of a unique axis given two pairs of points, as in step 3 above, for a 3D lattice.
Problem: Given two pairs of points on the unit sphere S 2 ⊂ R 3 , that is P 1 , P 2 and Q 1 , Q 2 .Assume we know there exists a rotation R with axis through the origin such that P 2 = R(P 1 ) and Q 2 = R(Q 1 ).Can we find the axis of this rotation by just knowing these four points?What about the angle of this rotation?Let us start with the following definition: Definition 4.1.A great circle of S 2 is the intersection of S 2 with a plane through the origin.For every great circle C ⊂ S 2 there are two well defined unit vectors v C , − v C which are perpendicular to the plane defining this great circle.There is no unique way which one of them to denote by v C and which one of them to denote by − v C , but we will use this notation since our arguments are essentially independent of this choice.
Note that any pair of two different great circles intersects in precisely two antipodal points.The following lemma states how we can find these two antipodal points: Lemma 4.2.Let C, C ⊂ S 2 be two different great circles with pairs of unit normal vectors ± v C and ± v C .Then the intersection points C ∩ C are precisely given via a normalised vector product by Proof.Let P, P ⊂ R 3 be planes through the origin such that C = P ∩S 2 and C = P ∩S 2 .Then ± v C and ± v C are unit normal vectors of P and P .Observe that any vector perpendicular to v C lies in P and any vector perpendicular to v C lies in P .Note that the vector product v C × v C is perpendicular to both v C and v C by the definition of the vector product.Moreover, we have v C × v C = 0 since the planes P and P are different and, therefore, their unit normal vectors v C and v C are linearly independent.Since v C × v C = 0 is perpendicular to both v C and v C , this vector must lie in both P and P , that is, in their intersection.Normalising this nonzero vector leads to an intersection point of the unit circles C and C in P and P .We have two intersection points and the other intersection point must be its antipodal point in S 2 .Now we state the first result: Theorem 4.3.Let P 1 , P 2 ∈ S 2 and Q 1 , Q 2 ∈ S 2 be two pairs of points such that there exists a rotation R centred at the origin with P 2 = R(P 1 ) and Q 2 = R(Q 1 ).Then the axis A of R is given as R a with a unit vector a ∈ S 2 ( a is determined only up to sign), and we can find a by the following normalised vector product: in the case that the vectors P 2 − P 1 and Q 2 − Q 1 are linearly independent.
Proof.The possible axes of rotations mapping P 1 into P 2 are determined by the great circle perpendicular to P 2 − P 1 .The possible axes of rotations mapping Q 1 into Q 2 are determined by the great circle perpendicular to Q 2 − Q 1 .If these two great circles are different, then they intersect in precisely two antipodal points, and these points must the unit vectors of the axis of the rotation R. The Lemma tells us that we can find these intersection points by computing .
We have a little problem in the case that the vectors P 2 − P 1 and Q 2 − Q 1 are parallel since, in this case, the potential axes mapping P 1 into P 2 are determined by the same great circle as the potential axes mapping Q 1 into Q 2 .The solution in this case is given by the following result: Theorem 4.4.Let P 1 , P 2 ∈ S 2 and Q 1 , Q 2 ∈ S 2 be two pairs of points such that there exists a rotation R centred at the origin with P 2 = R(P 1 ) and Q 2 = R(Q 1 ).Assume that the vectors P 2 −P 1 and Q 2 −Q 1 are parallel, so that we cannot apply the previous theorem.In this case, the axis A of R is given as R a with a unit vector a ∈ S 2 ( a is determined only up to sign), and we can find a by the following normalised vector product: Proof.If P 2 − P 1 and Q 2 − Q 1 are parallel, then the unique plane through the origin perpendicular to these two vectors is a reflection plane (mirror), mapping P 1 to P 2 and Then the unit vectors defining the axis of the rotation R are given as the intersection points of the great circles through P 1 , Q 1 and through P 2 , Q 2 .A normal vector to the great circle through P 1 , Q 1 is given by P 1 × Q 1 and a normal vector to the great circle through P 2 , Q 2 is given by P 2 × Q 2 .Now the results follows again via a direct application of the Lemma above.
Once we have the axis of rotation, it is not difficult to determine the angle of rotation for the pair P 1 , P 2 and for the pair Q 1 , Q 2 with respect to this axis and, since both pairs of points are associated to the same rotation, the angle of rotation must agree whether we compute it with the help of the pair P 1 , P 2 or with the help of the pair Q 1 , Q 2 .In fact we have the following: Theorem 4.5.Let P 1 , P 2 ∈ S 2 and R be a rotation centred at the origin with P 2 = R(P 1 ).Assume that a ∈ S 2 is a unit vector representing the axis of this rotation and that P 1 is not parallel to a. Then the angle θ of the rotation R is given by where "•" is the standard inner product in R 3 .
Proof.To find the angle θ of the rotation, we need to find the angle between the two planes E 1 and E 2 , where E i is the plane spanned by the vectors P i and a.This angle is the same as the angle between the unit normal vectors to the planes E 1 and E 2 .These unit normal vectors are given by a×Pi a×Pi .The result follows by the fact that a × P 1 = a × P 2 .

Mathematical deductions about 2D lattices allowing twinning
The algorithms described in Section 4 rely on finding lattice points whose distances from the origin are equal.In 2D, finding more than two points in the lattice at the same radius from the origin allows us to identify a rotation through an angle θ ∈ (0, π) such that some of the lattice away from the origin is mapped onto itself6 .A natural question to ask is: are there some lattices for which that are no non-trivial rotations which map one lattice point to another?In the context of crystallography, the question is whether there exist lattices that do not allow twinning to be detected from a single lattice.Without loss of generality, we consider the lattice Γ = Z a + Z b, generated by linearly independent basis vectors a = (1, 0) and b We denote this region in the upper half plane by M , that is (5.1) In this section we are concerned with determining which b ∈ M admit non-trivial pairs (m 1 , n 1 ), (m 2 , n 2 ) ∈ Z 2 , with (m 1 , n 1 ) = ±(m 2 , n 2 ), such that (5.2)

2D lattices allowing precise twinning
Rearranging (5.2), we find that the components of b must satisfy (5.4c) Note that α, β, γ ∈ Z are integers associated to a pair (m 1 , n 1 ) = ±(m 2 , n 2 ).The lattices We now distinguish two cases: (5.6) The condition γ = 0 means n 2 = ±n 1 .In the case (5.7) This case produces vertical rays at all rational points since we can choose any denominator n ∈ Z \{0} and and any numerator m 1 + m 2 ∈ Z we wish to have.Note that we only consider b Case 2: γ = 0: In this case, we can rewrite (5.3) as follows: which is the equation of a circle with centre (−β/γ, 0) and radius (5.10) We can rewrite this radius in terms of the pairs (m 1 , n 1 ) and (m 2 , n 2 ), and obtain the x-coordinate of the centre as (5.12) Note that both the centre of the circle of solutions b and its radius are rational numbers.
Let us understand the end-points of these semicircles hitting the x-axis perpendicularly.These end-points x 1 , x 2 ∈ R are given by x ± r and we obtain This shows that the end-points x 1 , x 2 of any of these circular semicircles meeting the x-axis perpendicularly are always rational.Interestingly, also the converse is true: each such semicircle with rational end-points appears.Namely, for any pair of rational points To see this we rewrite the ratios as 2a 2b and 2c 2d and we need to find This can be guaranteed by solving the linear equations The solution is given by We now return to the question we posed at the start of Section 5: are there some lattices for which that are no non-trivial rotations which map one lattice point to another?We have found the set, S, of all vectors b ∈ M which correspond to lattices that allow twinning: They are the union of all the vertical rays with rational x-coordinates with all semicircles with rational end-points on the x-axis, intersected with M .
The set of all vertical rays with rational x-coordinates is dense in R 2 , and correspondingly S is dense in M .We conclude that the set of lattices which allow twinning is dense in the set of all lattices in R 2 .Furthermore, note that the set S is the union of a countable number of curves.Each curve has Lebesgue measure zero, meaning S also has Lebesgue measure zero.We conclude that almost all lattices Γ ⊂ R 2 do not allow twinning.
This result is very interesting from a mathematical standpoint.In the following sections we will consider the necessary adjustments to make this result applicable to the experimental results from crystallography.

Truncated lattices
In the context of crystallography, the question we really want to answer is, for a set of basis vectors for a lattice can we determine a twin law using the algorithms from Section 4. So far we have not taken into account two important pieces of information about the experimental data: (1) Firstly, the intensity observed in real xrd data decays as the distance from the origin increases.This means that twinning can only reasonably be detected within some radius, N , of the origin.
(2) Secondly, the data measured in an xrd experiment is noisy, so it would be unrealistic to only seek points which are exactly the same distance from the origin.We should allow for some tolerance.
We start by addressing the first of these two points.Previously we found a set S ⊂ M of vectors b which satisfy (5.2).Each of the lattices, Γ, corresponding to each vector b ∈ S has the property that the lattice point m 1 a + n 1 b can be rotated non-trivially into a different lattice point m 2 a + n 2 b.However, for some lattices, these points could lie outside of the truncated lattice, (5.13) we are interested in.We therefore seek the set, S N , of basis vectors b ∈ M which satisfy For the following rough estimates of m 1 , n 2 , m 1 , n 2 we utilise the facts that a = (1, 0) and b is taken from the set M in (5. (5.15) for i ∈ {1, 2}.From this formula, it follows that |n i b 2 | ≤ N which implies Again using (5.15) we have Thus, we restrict the possible b which allow twin laws by reducing the rational values at which the above rays and circles can appear.Possible b which allow twin laws are contained in a finite number of semicircles and rays.These rays are restricted by b 1 = − m1±m2 2n and semicircles restricted to those semicircles with end points x 1 = m2−m1 n1−n2 and x 2 = − m1+m2 n1+n2 .The denominators of these are bounded by the bounds on n 1 , n 2 , so there is a finite number of such semicircles and rays.
We have found that S N of b-vectors corresponding to truncated lattices allowing twinning is not dense in M .This means that there is an increasing set of lattices where it is impossible to detect twinning within decreasing radii N of the origin.Knowing which lattices do not allow twinning means that time can be saved by not running a potentially computationally intensive algorithm to find a twin law that simply can not exist.

Including tolerance
So far in this section we have determined the lattices in 2D where it is possible to rotate a point onto another point exactly.However, experimental data is often noisy, and so it is unrealistic to expect data from an xrd to lie exactly the same distance from the origin.We therefore seek lattices where it is possible to rotate the lattice through an angle of θ ∈ (0, π) and have two points almost overlap.By this we mean there exists two points in the lattice, Γ = Z a + Z b, satisfying the relationship where > 0 is the experimental tolerance.As before (m 1 , n 1 ) = ±(m 2 , n 2 ) ∈ Z 2 .In Section 5.1, we found that the set S of b vectors which allow for exact overlaps is dense in M .We therefore expect that for each vector b ∈ M , we will be able to find two points satisfying (5.16).We denote the minimum distance from the origin of a point satisfying (5.16) by R( b).We discretize the space M , and in Figure 7 we plot the calculated R( b) values, using a lighter blue as R( b) increases.As we have assumed without loss of generality that b ≥ 1, the data within the unit circle should be ignored.
In Figure 8, we plot those b vectors which lead to almost exact overlaps within the truncated lattice Γ N , with the coordinates of the b vectors with R( b) ≤ N shown in blue, and R( b) > N shown in yellow.We recognise the line segments and circular arcs predicted in Section 5.2.We are able to repeat this calculation for different values of N and , and the resulting data could be used as a look up table to see if it is possible to detect twinning in an experimentally determined lattice.The numerical investigation in this section in particular is amenable to extension into 3D.

An interesting connection with hyperbolic geometry
There is a stunning agreement between our Figure 8a and a figure from [3], see Figure 9.This is, indeed, not a coincidence.The curves in that paper represent closed geodesics in the modular surface and the curves appearing in our figure can be understood as geodesics starting and ending in cusps of the modular surface.The modular surface is well studied and an object in a very different geometry, called Hyperbolic Geometry, Let us provide a brief overview over hyperbolic geometry: Our model space is the upper half plane Here it is better to work in the complex plane C instead of R 2 .In this world, smooth curves α : [a, b] → H 2 have a different length from their usual Euclidean length, and we refer to this length as their hyperbolic length: where α(t) = x(t)+iy(t).This new length definition allows us to introduce a new distance function between points of H 2 , that is, a function d : There is even an explicit formula for this hyperbolic distance function, namely makes H 2 a metric space.This metric space is called the hyperbolic plane or, since there are also other models of this geometry, the upper half space model of the hyperbolic plane.It is a space of constant negative Gauss curvature minus one.It is the counterpart to the unit sphere which is a space of constant positive Gauss curvature plus one.Let us briefly discuss the notion of a geodesic: A geodesic in H 2 is a curve α : I → H 2 (with I ⊂ R an interval) of shortest between any two of its points.In this respect, a geodesic is a straightforward generalisation of a straight Euclidean line segment in the Euclidean plane R 2 .In other words, if t 1 , t 2 ∈ I with t 1 < t 2 and if β : [t 1 , t 2 ] → H 2 is the restriction of α to this subinterval, then we have it turns out that all geodesics of H 2 have a very specific shape: they are either vertical Euclidean line segments in H 2 (that is, subsets of the form {x + iy : y ∈ [y 1 , y 2 ]}) or they are segments of Euclidean circles which meet the x-axis perpendicularly.Here it becomes apparent that the curves appearing in our pictures may be geodesics.The group acts via Möbius transformations on H 2 , that is, These Möbius transformations are isometries of the hyperbolic plane H 2 , that is, we have for z, w ∈ H 2 and A ∈ SL 2 (R): There is a natural subgroup of SL 2 (R), namely the group and the role of M for the SL 2 (Z)-action on H 2 is the same as the role of Z for the Γ-action on R 2 .)Note that SL 2 (Z) is a discrete subgroup of SL 2 (R), the latter acting by orientation preserving isometries on H 2 .Such discrete subgroups are generally referred to as Fuchsian groups.
The quotient H 2 /SL 2 (Z) is called the modular surface and it is constructed from the fundamental domain M by specific side identifications, namely, by identifying the points of the side −1/2 + iy with the corresponding points of the side 1/2 + iy and also by identifying the point −x + iy with −1/2 ≤ x ≤ 0 and x 2 + y 2 = 1 with the point x + iy.(This is the same process as the construction of the torus R 2 /Γ from the unit cell Z: There we identify the boundary points v b with a + v b and the boundary points u a with u a + b.)It turns out that the modular surface H 2 /SL 2 (Z) is a noncompact space with one cusp corresponding to the point infinity of the hyperbolic plane H 2 .Geodesics in the modular surface are just projections of geodesics of the hyperbolic plane into this modular surface and geodesics which both start and end in the cusp of the modular surface are precisely projections of geodesics in the hyperbolic plane H 2 which are either Euclidean semicircles starting and ending at rational points of the x-axis or straight vertical Euclidean rays, starting at rational points of the x-axis.This completes the explanation of this surprising connection between two different geometries.For more information about hyperbolic geometry we refer readers to [1] and Fuchsian groups of which SL ( 2, Z) is a specific example to [4].
Isn't this a nice ending of our mathematical journey, which started by an innocent twinning question?

Conclusion
Throughout this work, we have investigated various paths towards a better twin finding algorithm.In Section 2, we discussed how fake data would be generated to closely mirror true XRD data.In Section 4, we discussed an algorithm based on the rules for rotations outlined in 3. Finally, Section 5 investigated the possibility of 2D lattices which do not allow any possibility of twinning, and found a curious link with hyperbolic geodesics.
We are confident that the results presented here represent a significant development for the post-hoc detection of twin laws, and are excited to see what further developments can be made.

Outlook
There is much within this work which can be expanded upon.The current search for 3D overlaps of rotations in Section 4 is able to present a set of ellipses, of which points where many pass close together will give a very viable twin law, but we have no present algorithm for simply finding such points, and rely on a crude search of space.If intersections (and particularly almost-intersections) of such ellipses could be found in a more analytical fashion, we may gain a dramatic time improvement in the detection of twin laws.
Section 5 considered lattices for which twinning would not be detectable in 2D.If this were extended to 3D space, it could provide an incredibly simple check to prevent testing for twins when there are none present.Additionally, an expansion of this test to look for not only a single overlap within the threshold, but for a particular number, could also encourage the 'ruling out' of twin laws, as the detection of a twin law is typically the result of many such overlaps.

Figure 2 .
Figure 2. A pyrite crystal with three components

Figure 3 .
Figure 3. Two-dimensional fake data for a square lattice.The size of the points are proportional to the simulated intensities.(a): Intensities for no-twinning.(b): The same case as (a), but with a twinned component with twinning angle 36.9°andβ = 0.1.(c): The lattice rotated by the twinning angle, with the green dots showing the residual intensities between the twinned and no-twinning casces.
(a) Azimuthal angle = 30°( b) Azimuthal angle = 45°F igure 4. Three-dimensional fake data for a cubic lattice, with all points the same size for visual clarity.(a): The data at an azimuthal angle of 30°, where the data appears to be unstructured.(b): The data at an azimuthal angle of 45°, where the some of the lattice structure can be seen in the vertical lines Figure 5. Axes of rotation in 3D

Figure 6 .
Figure 6.An example of the output of Step 4 for a 2D lattice.In this case, one of the angles revealed by the Step 4 process is 67.38 • .The underlying lattice is shown in blue, and the lattice points that are rotated onto each other by this rotation are shown in red.
a + Z b corresponding to the solutions b = (b 1 , b 2 ) of (5.3) have the property that the lattice point m 1 a + n 1 b can be rotated into the lattice point m 2 a + n 2 b.

Figure 7 .
Figure 7. Minimum radius, R( b), at which a point in the lattice Γ = Z a + Z b, a = (1, 0), b ≥ 1 satisfies (5.16).Lighter shades of blue denote large values of R( b).Inside the unit circle should be ignored.The tolerance is set to = 0.005.

Figure 8 .
Figure 8.The coordinates of b ∈ M with R( b) ≤ N are shown in blue.Inside the unit circle should be ignored.The tolerance is set to = 0.005.
Figure 9.The distribution of G 377 projected on the fundamental domain of SL 2 (Z)/H from [3, Figure 2]