Catastrophe theory classification of Fermi surface topological transitions in two dimensions

We classify all possible singularities in the electronic dispersion of two-dimensional systems that occur when the Fermi surface changes topology, using catastrophe theory. For systems with up to seven control parameters (i.e., pressure, strain, bias voltage, etc), the theory guarantees that the singularity belongs to to one of seventeen standard types. We show that at each of these singularities the density of states diverges as a power law, with a universal exponent characteristic of the particular catastrophe, and we provide its universal ratio of amplitudes of the prefactors of energies above and below the singularity. We further show that crystal symmetry restricts which types of catastrophes can occur at the points of high symmetry in the Brillouin zone. For each of the seventeen wallpaper groups in two-dimensions, we list which catastrophes are possible at each high symmetry point.


I. INTRODUCTION
Valence electrons in crystalline solids are described by Bloch states with a dispersion relation n (k) between energy and crystal momentum k, with n denoting a set of discrete indices such as band, spin, etc. Within this construction, the density of states (DOS) as a function of energy or momenta plays an important role in the calculation of physical properties such as heat capacity and magnetic susceptibility. While this is effectively a non-interacting picture, the electronic dispersion serves as an input for the treatment of interactions. In particular, an enhancement of the DOS at the Fermi level can strengthen the role of electron-electron interactions and lead to electronic instabilities.
In one and two dimensions, a divergence of DOS can occur due to the presence of one or more critical points in k-space where ∇ (k) = 0 [1] (from here on we drop the index n). The simplest such singularity in two dimensions, the van Hove singularity (vHs), occurs due to regular saddle points and causes a logarithmic divergence of the DOS. The Hessian at a regular saddle point is non-degenerate [i.e., det(∂ i ∂ j ) = 0] and the dispersion is quadratic to the lowest order [ (k) ∼ k 2 x − k 2 y ]. However there also exist higher order singularities which correspond to power law divergence of the DOS. These occur due to higher order critical points at which the Hessian is degenerate, i.e. non-invertible, and the dispersion needs to be expanded beyond quadratic order.
Singularities and the associated divergence of DOS are a signature of Fermi surface topological transitions [2,3] in which the Fermi surface undergoes a change in topology from electron type to hole type across the critical energy, with two or more branches of the Fermi surface touching at the critical point in a singular way. Historically, when Lifshitz [4] first studied the change of Fermi surface topology, he dealt with two cases: the appearance or collapsing of a neck and the appearance or collapsing of a pocket in the Fermi surface. The neck-collapsing case is the ordinary vHs, with the Fermi surface locally consisting of a pair of intersecting straight lines. These two types of Fermi surface topological transitions have been observed, along with their non-trivial consequences due to the presence of interactions, in a wide range of quantum materials including cuprates, iron arsenic and ferromagnetic superconductors, cobaltates, Sr 2 RuO 4 , heavy fermions [5][6][7][8][9][10][11][12][13].
Higher order singularities display more exotic Fermi surface topological transitions. They have recently been associated with phenomena such as the unusual Landau level structure and tripling of de Haas−van Alphen and Shubnikov−de Haas oscillation periods in biased bilayer graphene [2], the non-trivial thermodynamic and transport properties in Sr 3 Ru 2 O 7 [3], correlated electron phenomena in twisted bilayer graphene near half filling [14] and the so called supermetal with diverging susceptibilities in the absence of long range order [15]. In the present work we develop a classification scheme for Fermi surface topological transitions and their associated DOS divergence using catastrophe theory. That catastrophe theory is an appropriate framework for examining higher order critical points in electronic bands as was first suggested by [2]. It has been also applied to other branches of physics (e.g. Ref (16), including optics [17] and molecular physics [18]).
Catastrophe theory deals with real valued functions of open subsets of R n , where R is the space of real numbers, and makes a distinction between state variables and control parameters. In the context of electronic bands the function of interest is the electronic dispersion and we identify the components of the crystal momentum as the state variables and hopping strengths, chemical potential, etc as control parameters which can be tuned externally, for example, by applying pressure, strain, bias voltage, etc. At a critical point, the gradient of the function with respect to the state variables vanishes. Since we focus on two dimensional systems, we are restricted to exactly two state variables, k x and k y . With two state variables, catastrophe theory guarantees that a system with seven or less control parameters is typically equivalent to one of seventeen standard types of catastrophes, each of which corresponds to a unique higher order singularity. The higher order singularities are indexed by three positive integers: the corank, determinacy and codimension. The classification by these numbers is unique except for certain arXiv:1910.10183v1 [cond-mat.str-el] 22 Oct 2019 degenerate cases, which we show that in two-dimensions can be further distinguished by the winding, i.e., the number of times the electronic dispersion changes sign along a contour encircling the critical point.
By tuning the control parameters, one can reach the higher order critical points corresponding to different catastrophes. In [3], the fourth order saddle was identified with the unimodal parabolic singularity X 9 while in [2] the monkey saddle was identified with the elliptic umbilic catastrophe. In the same spirit, we identify the higher order singularity treated in [14] and [15] with the cusp catastrophe.
In the example in [3], the elliptic umbilic catastrophe occurs on a point with 2π/3 rotational symmetry, while in that of [15] the cusp occurs on a point with π rotational symmetry. As a general principle, symmetry constrains a point to be critical and serve as the center where ordinary critical points merge. Only catastrophes consistent with the symmetry can occur at such a point. Another feature of high symmetry points in the Brillouin zone is that they can host otherwise atypical higher order singularities that are not part of the standard seventeen. These facts call for further classification of the catastrophes that can occur at the high symmetry points. We present such a classification for the Brillouin zones corresponding to the seventeen wallpaper groups (no relation to the seventeen catastrophes).
The paper is organized as follows: In Sec. II we introduce the language of catastrophe theory through a simple example of a tight-binding model that displays a higher order singularity. This is followed by the classification of singularities in electronic bands in Sec. III. The stability of this classification to small perturbations is discussed in Sec. IV. In Sec. V, we explain the connection between high symmetry points in the Brillouin zone and higher order singularities. We then present a classification of the singularities allowed at the high symmetry points in the Brillouin zones corresponding to the seventeen wallpaper groups. In Sec. VI, we give a practical method for determining the type of singularity given its Taylor expansion and illustrate the method with a sample calculation. Finally we conclude the discussion in Sec. VII by briefly sumarizing the scope of the work and setting the context for future work on the treatment of line singularities and the effects of interaction.

II. A SIMPLE EXAMPLE
In this section we introduce the languge of catastrophe theory through a simple tight binding model that displays a higher order singularity.
The lattice pertaining to the model is layered with two sublattices A and B, colored respectively by black and grey in Fig. 1. In theŷ direction, we have only AA and BB nearest neighbor (NN) hoppings of strength t 2 . In thex direction, we have A → B NN hopping of strength t 1 and imaginary A → A and B → B next nearest neighbor (NNN) hopping it . The imaginary hopping can be interpreted as arising either from spin-orbit coupling or a staggered magnetic flux in theŷ direction.
FIG. 1. The lattice for the Hamiltonian treated in Sec II is shown above. It is a square lattice with a two site basis of atoms colored as black (A) and grey (B). The AB nearest neighbor (NN) hopping in thex direction has a strength of t 1 while the AA and BB NN hoppings in theŷ direction have strength t 2 . There is also complex next nearest neighbor AA and BB hopping it (with direction as depicted in the figure). These hopping strengths can be tuned so that the energy bands of this Hamiltonian yield the simplest of the higher order singularities, the fold.
Choosing a unit cell containing one A atom and one B atom as shown in Fig. 1, the Hamiltonian is given by: where c A/B,R annihilates an A/B type Fermion in the unit cell located at R. Diagonalizing this Hamiltonian in k-space yields two bands: ± (k) = ±2t 1 cos (k x /2) − 2t sin (k x ) + 2t 2 cos (k y ). We now expand + (k) to cubic order around the point (π, 0) in the first Brillouin Zone (BZ): where α = 2t − t 1 and β 3 = t 1 /4 − 2t are independent. To obtain a critical point, we need to remove the k x term by tuning α to zero. We assume β = 0 and t 2 = 0. This guarantees that we can always rescale (k x , k y ) → (k x /β, k y / √ t 2 ) and write t = α/β to obtain: Thus, up to a rescaling of coordinates and a constant energy shift (2t 2 ), we effectively have just one control parameter t in y + tk x for three cases: (a) t < 0, (b) t = 0 and (c) t > 0. When t < 0 there is an ordinary maximum and an ordinary minimum. As t is reduced to 0, these merge at the origin to give a higher order singularity. For t > 0 there is no longer a critical point. this dispersion.
In Fig. 2 we illustrate the behaviour of the energy contours and the Fermi surface for three distinct cases. When t < 0, the Fermi sea is topologically non-trivially connected, with an island of positive energies (small red region) in the middle of the Fermi sea (blue region). There are two ordinary critical points in k-space, demarked with thick dots: a maximum in the dispersion (the summit of the red island) and minimum (bottom of the blue Fermi sea). Upon increasing t to 0, a higher order critical point (k 3 x − k 2 y ) appears suddenly. This higher order critical point at t = 0 corresponds to the merger of the two ordinary critical points seen for t < 0, as well as to the Fermi surface topological transition in which the island disappears and the Fermi sea becomes trivially connected. When t > 0, the dispersion is no longer singular, and the single component Fermi sea persists.
At the higher order critical point (t = 0), the DOS diverges as | − 2t 2 | −1/6 when → 2t 2 . By including higher order terms in the dispersion we still obtain the same power law divergence to the lowest order and the Fermi surface topology does not change. The value of the exponent is a property of the type of catastrophe, and cannot be changed by smooth coordinate transformations, as we show in Appendix A. For the singularity in this simple example, the cubic truncation of the Taylor series is sufficient.
As the above example shows, we can often tune some control parameters in a system to obtain a higher order critical point. By a suitable smooth transformation of the coordinates and the control parameters, we can put the dispersion in a convenient form that contains the higher order singularity, to which polynomial terms modulated by the control parameters are added (Eq. 3). For a given higher order singularity, there is usually a most general such expression, known as the universal unfolding, containing only a finite number of effective control parameters modulating polynomial perturbations (t k x in the example). These polynomial perturbations represent directions in polynomial space which are not tangent to the orbit of the higher order singularity under smooth coordinate transformations. The number of such missing directions (or equivalently effective control parameters) is the codimension of the singularity. Finally, these higher order singularities are well described by truncating the Taylor expansion to a certain lowest order, known as the determinacy. In the example above, the determinancy is 3, justifying the truncation to cubic order.
In addition to the codimension and determinacy, catastrophy theory utilizes another index, the corank, which is the number of zero eigenvalues of the Hessian matrix of the function at the critical point. (See Appendix C for a self-contained review of catastrophe theory). While these three indices can classify singularities in any dimension, they are not sufficient to resolve the degeneracies. In the particular case of functions of two variables, (for instance k x and k y in the case of dispersions of two-dimensional electronic systems), we find that these degeneracies can be resolved by introducing a fourth index, the winding number, which counts the number of times the function changes sign in a small loop around the origin.
In terms of these four indices (codimension, determinacy, corank, and winding), we can enumerate all possible singularities that the dispersion of a two-dimensional electronic system would typically present. These singularities, along with the associated power laws in the DOS and universal ratios of their prefactors, are listed in the next section.

III. CLASSIFICATION OF HIGHER ORDER SINGULARITIES IN ELECTRONIC SYSTEMS
Using catastrophe theory we can classify all possible higher order singularities that can occur in the electronic dispersion of two-dimensional electrons. If there are seven or less control parameters in the system, catastrophe theory guarantees that only higher order singularities with codimension (or cod) 7 are typically likely to occur. Typicality here has a simple but precise mathematical meaning, that is best illustrated by an example. Typically a set of three equations in two variables is over determined and does not have a solution (two straight lines in three dimensional space do not typically intersect), while a system of two equations in two variables typically has a unique solution, at least locally (two straight lines in two dimension typically intersect at a point).
There are 17 singularities with cod 7, which we list in Tables I, II and III. We also include one other atypical singularity, the X 9 , which has codimension 8. It would normally require 8 control parameters, but crystal symmetry ensures that some of the constraints are automatically satisfied. In fact it is the lowest codimension singularity that respects four-fold rotational symmetry. The tables are organized by the four integer indices (cor, cod, det, w) corresponding to corank, codimension, determinacy, and winding.
As mentioned earlier, at an ordinary critical point in two dimensions, the DOS has a logarithmic divergence. At higher order critical points, the DOS diverges in a power law fashion, often with different coefficients D + and D − as we approach the critical energy from above and below: In Appendix A we show that both the exponent and the ratio of the coefficients are universal in that they are preserved under smooth coordinate transformations. The power law dependence for dispersions with two monomial terms can easily be extracted by scaling k x and k y appropriately, as mentioned in [15]. But this procedure by itself does not give the prefactors and one would still have to convert the DOS integral into an integral over constant energy contours to evaluate D ± . This is also done in Appendix A. In the tables, we list for each catastrophe both the exponent γ controlling the divergence of the DOS at the critical point and the universal ratio D + /D − .

IV. ROBUSTNESS OF THE CLASSIFICATION
Given that higher order critical points are somewhat rare and may require tuning of parameters to obtain, it is natural to ask how robust the classification scheme presented above is. This is particularly relevant in real systems which are prone to imperfections manifesting as small perturbations, motivating the question: can small perturbations change the type of higher order singularity that can occur on tuning parameters? The answer is no. This is because the universal unfolding contains all possible perturbations to the higher order singularity. For any sufficiently small perturbation, we get the same catastrophe, up to a smooth coordinate transformation. Only the critical value of the tuning parameters for which we obtain the higher order singularity, gets shifted. This property is referred to as stability in the literature. It is important to note that stability is guaranteed only for singularities with codimension 7, viz the seventeen listed above. For codimension 8 and higher we encounter the so called unimodal singularities which possess a continuous parameter in their dispersion. Two singularities with different values of the continuous parameter, however close by, will in general not be smoothly equivalent so that a small perturbation can take us to an altogether different singularity. But as mentioned earlier these are atypical in systems where seven or less parameters are tuned.

V. HIGH SYMMETRY POINTS AND CATASTROPHES
As seen in Sections II and III, within a Bloch-band description, a point in k-space can morph into a higher order critical point when the parameters in the Hamiltonian are tuned. While this is may be hard to achieve at an arbitrary point, high symmetry points in the Brillouin zone are particularly easy to tune into higher order critical points. Recall that for a point to be critical, we need both components of the Jacobian at the point to vanish. For higher order critical points, we further need the determinant of the Hessian to vanish (or equivalently at least one of the eigenvalues needs to be zero). If a reflection axis passes through the point, the component of the Jacobian perpendicular to the axis is already constrained to be zero, yielding a simplification. Alternatively, if the point is a center of a non-trivial rotation, the entire Jacobian will vanish so that we have to tune only the determinant of the Hessian to zero (see Appendix B for details). Thus it is relatively easier to obtain higher order singularities at points of high symmetry by tuning parameters in the Hamiltonian.
The Brillouin zone of the wallpaper group p4 shown along with the high symmetry points. The point group of p4 is generated by π/2-rotation which gets carried over to k-space as a symmetry that acts about the origin. By combining the π/2 and π rotations with reciprocal lattice translations we can identify the other high symmetry points. At each of the high symmetry points, only dispersions which respect the symmetry are allowed. This places a huge restriction on the higher order singularities that are likely to occur. In the insets we identify the catastrophes likely to occur at each of the high symmetry points.
Having stressed the importance of high symmetry points in k-space, we remark that they can be identified in a systematic way (see for instance [19]). Given a lattice with a certain symmetry group, the dispersing bands have the symmetry of the point group. If the Hamiltonian is time reversal    Table I they have saddle like nature with a non-zero winding number number that counts the number of times the dispersion changes sign along a closed contour enclosing the origin. However, unlike the cuspoids wherein those with even degree have π rotation symmetry, most of the umbilics do not possess any rotational symmetry. The elliptic umbilic (monkey saddle) is an exception in that it has 2π/3 rotation symmetry.
The original Thom's theorem dealt with seven catastrophes. It was first extended to twelve and then to seventeen catastrophes. The last five catastrophes are documented above along with the X 9 singularity which is the singularity of lowest codimension that is allowed by π/2 rotation. Unlike the rest, it is unimodal: different values of the parameter c in the dispersion correspond in general to inequivalent singularities with different windings. The Fermi surface shown above is for c = −3, a generalization of the monkey saddle. invariant (TRI) and there is no spin orbit coupling, we have an additional inversion symmetry in k-space. By combining the point group symmetries with reciprocal lattice translations we can find the high symmetry points in the first BZ. We illustrate this for the wallpaper group p4 in Fig 3. The point group of p4 is C 4 , which is generated by π/2-rotation. This point group gets carried to k-space wherein it acts about the Γ point (the origin in k-space). Γ is thus a center of four fold rotation. We then combine the π/2-rotation and π-rotation with various reciprocal lattice translations.
Combining rotations with translations simply shits the center of rotation. This procedure then yields the M and X points as centers of 4-fold and 2-fold rotations respectively in the first BZ.
At the high symmetry points, only singularities consistent with the symmetry of the point are allowed. For example, at the Γ and M points in the above example, only singularities which are invariant under π/2 rotation are allowed to occur. X 9 is the singularity with lowest codimension that satisfies this property and is thus the one that is typically likely to occur. At the X point, only the cuspoids A 2n−1 with the standard form k 2n x ± k 2 y for n = 1, 2, 3 are allowed to occur. This is because these are the only ones among the seventeen which are invariant under π rotation. (It can be shown that no smooth coordinate transformation can make the rest among the seventeen consistent with π rotation). In Table IV we apply this scheme for the Brillouin zones corresponding to the seventeen wallpaper groups in two dimension. (See Appendix F for a quick review of crystal symmetries and their consequences for bands).

VI. IDENTIFYING A SINGULARITY IN PRACTICE
In practice, the singularities need not occur in their standard form. Each of the 17 singularities forms a distinct equivalence class. (Equivalence here is up to a smooth, local coordinate transformation with a smooth inverse). So we generally expect some smooth equivalent of the standard forms to occur in the context of band dispersions. But it might be hard to identify the singularity from such a non-standard expression. While we could try exhibiting a smooth coordinate transformation that puts the singularity in its standard form, it is very hard to do that for a generic dispersion. This apparent complication is easily resolved by exploiting the fact that codimension, determinacy, corank and winding are invariant under smooth coordinate transformations. By computing these numbers for the given singular dispersion we unambiguously identify its type. Below we give an explicit method for computing determinacy and codimension.

A. Algorithm for computing determinacy
We now describe the method for computing determinacy and codimension and illustrate it with an example. The example has been chosen to illustrate the subtle aspects of Taylor expansions and how truncating the Taylor series to an incorrect order can yield wrong results which might appear sensible. We start by giving a few definitions and notations. The k th degree Taylor polynomial of f denoted by j k [f ] or simply j k f , is the truncation of the Taylor expansion of f at the origin to degree k. A function f is said to be k-determinate if any function having the same k th degree Taylor polynomial as f is smoothly equivalent to it. In particular of course, j k f is also smoothly equivalent to f . Determinacy is then the lowest k for which f is k-determinate.
The method given below is a simplification of the method given in [20]. It is also similar to the methods given in [21]. (See Appendix E 1 for an explanation of the method).
1. Choose a particular k to start. Determine the k th degree Taylor polynomial j k f .
2. Compute all polynomials of the form k m x k n y ∂j k f ∂k x and k m x k n y ∂j k f ∂k y with 1 m + n k and truncate them to degree k. There is only a finite number of such polynomials and we denote them by p 1 , · · · , p N . These can be thought of as vectors in the ploynimial space spanned by the set {k m x k n y ; 1 m + n k}.

Check if the system of linear equations
has a solution for each j = 0, · · · , k.
4. If a solution exists for each of the systems, f is kdeterminate. We keep reducing k by 1, repeating the algorithm for each k until we find the smallest k for which the systems of equations all have a solution. Call it k low and terminate the algorithm.
5. If solution does not exist for each of the systems, we keep increasing k by 1, applying the algorithm for each k until we find the first k for which all of the systems of equations have a solution. Call it k low and terminate the algorithm.
Once the algorithm terminates successfully and k low is determined, it can be shown that the determinacy of f is either k low or k low − 1. While this ambiguity cannot be resolved easily, it is not a problem as long as we determine the other indices unambinguously. This is because no two of the seventeen singularities whose determinacy differs by 1 have the same corank, codimension and winding.
Let us now give a method for finding the codimension. For k = k low , construct a matrix M by listing each of p 1 , · · · , p N , ∂j k f ∂k x , ∂j k f ∂k y as row vectors in the basis of the monomials {k m x k n y ; 1 m + n k}. Let rank(M ) denote the rank of M . The codimension is then (

2-fold rotation center (π)
3-fold rotation center (2π/3) 4-fold rotation center (π/2) 6-fold rotation center (π/6) reflection axis TABLE IV. In two dimensions, there are seventeen wallpaper groups which serve as symmetry groups of crystals. As a consequence of Bloch's theorem, the point group of the wallpaper group (mentioned above within parenthesis) acts about the origin in k-space. By combining the point group elements with the reciprocal lattice translations, we find the high symmetry points in the Brillouin zone and list the catastrophes that are allowed at each of these points. If the Hamiltonian has time reversal invariance and no spin orbit coupling, there is an additional inversion symmetry in k-space. Such cases are denoted above with a * (For example p1*).

B. A sample computation
Consider two functions f (k x , k y ) = −k 2 y + 6k 2 x k y − 8k 4 x and g(k x , k y ) = −k 2 y + 6k 2 x k y − 9k 4 x + k 6 x . On cubic truncation, both these functions yield h(x) = −k 2 y + 6k 2 x k y . This function arises for example in the context of the continuum model for twisted bilayer graphene in [14]. The Fermi surfaces of the three functions are depicted in Fig 4. Visually, they appear to share the same topology for the energy contours. In particular, the zero energy contour seems to consist of two curves intersecting tangentially at the origin. Naively this might motivate truncation to cubic order for computing the DOS among other things. We now demonstrate that this would lead to an incorrect result, at least for g, by showing that while f and h are equivalent to the cusp, g is equivalent to butterfly.
We first compute the determinacy of f . Let q n denote the truncation of the polynomial q to degree n. Choosing k = 4 we list some of the polynomials p i for this function: The common factors will not affect the results and we have dropped them. It immediately follows that 1 12 This verifies that f is equivalent to cusp. Proceeding in the same way we can show that det[g] = 6 rather than 4 and that it is equivalent to butterfly. This clearly shows that f and g are different singularities inspite of the similarity of their Fermi contours.
Finally it can be shown that h is also equivalent to the cusp. This naturally leads to the question as to why g is also not equivalent to the cusp given that its cubic truncation is h, which is equivalent to the cusp. The answer lies in the definition of determinacy. The relation det[h] = 4 only guarantees that functions whose fourth degree Taylor polynomials match that of h are equivalent to it and therefore the cusp. This importantly requires such functions to have no quartic terms in their Taylor expansion as h does not. Since g fails to satisfy this property, it is not guaranteed to be equivalent to h. This also means that from the point of view of the mathematical theory, the truncation of f to cubic order is also not justified although in this specific instance it worked.

VII. SUMMARY AND OUTLOOK
In this paper we presented a classification of possible Fermi surface topological transitions in two spatial dimensions that take place via higher order singularities. At these higher order singularities, the DOS diverges as a power law, in contrast to the weaker logarithmic divergence at van Hove singularities. When the number of control parameters is less than or equal to seven, we can list all the seventeen possibilities classified by catastrophe theory. The critical exponents controlling the divergence of the DOS, as well as the ratio of the prefactor or amplitudes for energies on both sides of the transition, are universal properties associated to the specific type of catastrophe.
Smooth coordinate transformations cannot change these exponents or ratios, as well as the topology of the critical Fermi surface right at the singularity.
We further showed that by tuning the parameters in the Hamiltonian, we can obtain higher order singularities. They are especially easy to achieve at high symmetry points in the Brillouin zone. At the same time, the high symmetry points also restrict the types of singularities: only those consistent with the symmetry can occur. Motivated by this, we classified the singularities that are likely to occur at high symmetry points in the Brillouin zones corresponding to the seventeen wallpaper groups in two dimension.
These critical points of the electronic dispersion can serve as the starting point for the treatment of interactions, and can potentially lead to novel electronic phenomena. Besides simply facilitating ordering transitions because of the enhancement of the DOS, there are several more exotic possibilities −1  x k y shown within the range −1 < k x < 1, −1 < k y < 2. While the energy contours of f and g appear to have the same topology, they are actually not related by smooth deformation. The former corresponds to two parabolas that intersect tangentially and the latter consists of two cubic curves that intersect tangentially at (0, 0). The DOS diverges differently for both these dispersions although they give the same cubic truncation h. This in spite of the fact that h appears to produce correct topology of Fermi surface (two curves tangentially touching at an isolated higher order critical point) and a seemingly sensible power law divergence of DOS. This is shown to be the consequence of the fact that f and g are inequivalent and have different allowed order of truncation as indicated by their determinacy. that demand exploration. While the addition of weak interactions on smooth Fermi surfaces leads to a Landau Fermi liquid, it is not clear that singularities in the curvature of the Fermi surface and their associated power law enhancement of the electronic DOS may not lead to novel quantum critical behavior. In the singularities that we discussed, the singular point is not isolated (i.e., a point-like Fermi surface); there is a finite Fermi surface containing this point. So the analysis of interactions for certain singularities may require dealing with singular and non-singular parts of the Fermi surface. For instance, the singular points with their divergent DOS may enhance scattering and reduce quasi-particle lifetime; at the same time, the Fermi velocities vanish at those points. Therefore, properties like transport, which depend both on the Fermi velocities and quasiparticle lifetimes, may depend on the joint effects of the singular and non-singular parts of the Fermi surface. Understanding the effects of interactions on the higher order singularities is an interesting and largely open problem; while some progress has been recently made [2,3,14], much remains to be understood.
In addition to the understanding of the classes of Fermi surface topological transitions, it should be possible to apply the machinery of singularity and catastrophe theory to other aspects of electronic dispersions. Here we focused on point singularities, but it may be fruitful to develop a catalog of line singularities as well. For example, in tight-binding Hamiltonians for the diamond lattice, there are cases where the Fermi surfaces collapse to lines. A step towards these generalizations would be to start with polynomials in k x , k y that multiply Pauli matrices, instead of just scalars. The matrix structure should bring about square root singularities that are not present in the simple scalar dispersions. We believe that the study here presented is one step in understanding singular behavior in electronic systems by deploying the tools of catastrophe theory.
Note added: While we were completing this manuscript, we became aware of similar work also being completed by Noah F. Q. Yuan and Liang Fu. We point the reader to that reference as well. We thank Zhi-Cheng Yang for recognizing that the works were similar, and bringing the two groups into contact. The DOS in the n th band ignoring spin degeneracy is defined as [22] (A1) Where S n ( ) is the constant energy contour in the n th band. If for a particular energy, the constant energy contour contains a critical point, the integrand diverges at that point.
Let us drop n so that the given dispersion isẼ(k). Assume that under a suitable local smmothly reversible coordinate transformationk → k, the dispersion transforms to E(k), one of the standard types. The Jacobian determinant for the coordinate transformation can be Taylor expanded as: where J 0 = 0 since that condition is equivalent to the coordinate transformation being a local diffeomorphism. Since the required coordinate transformation will in general be local, we restrict the DOS integral to a bounded neighborhood D of the origin. This procedure will not affect the calculation of the divergent part of the integral since it arises due to the critical point at the origin. (In fact we can also extend the integral to the entire two-dimensional plane). The DOS now reads: For dispersions which are either homogeneous or take the form E(k) = a k m 1 x k n 1 y + b k m 2 x k n 2 y for non-negative integers m i and n i , it might be possible to scale k x → | | α k x and k y → | | β k y for appropriate positive rational numbers α and β so that For instance, for a homogeneous dispersion of degree n, α = β = 1/n while for E(k) = a k m 1 x k n 1 y +b k m 2 x k n 2 y , we set α = (n 2 −n 1 )/(m 1 n 2 −m 2 n 1 ) and β = (m 1 −m 2 )/(m 1 n 2 −m 2 n 1 ). (For the standard types which take the latter form, m 1 n 2 −m 2 n 1 = 0 so that this scaling is allowed). Performing this scaling in Eq. A3, we get: Since α, β, γ > 0 and n m, we have mα + (n − m)β > 0 so that in the limit → 0 the leading order divergent part of the integral is given by: where once again we choose −1 if > 0, +1 is < 0. We have extended the integral to the entire plane since this only adds only a finite error that does not change the coefficient of the leading order divergent part. Thus, we infer that the exponent γ and the ratio D + /D − are invariant under a smoothly reversible coordinate transformation. We compute D + /D − for the various catastrophes in the forthcoming sections.

Cuspoid catastrophes
The cuspoid catastrophes include the fold, cusp, swallowtail, butterfly, wigwam and star whose dispersions take the form: (k) = k n x − k 2 y . For the sake of notational simplicity, in what follows we replace k x → x and k y → y. We now derive the DOS for these.
Consider (x, y) = x n − y 2 . This gives y = ± √ x n − . Furthermore: . This gives ∇ = n 2 x 2n−2 + 4y 2 = n 2 x 2n−2 + 4x n − 4 , where we have used the dispersion to substitute for y 2 . Now let n be even. For < 0, the constant energy curves have vertices (0, ± | |) disperse upward/downward in the xy-plane. The DOS integral then becomes: x n − (A7) Where we have substituted for dl and ∇ in the integrand and simplified the expression. We have also accounted for the fact that the upward and downward dispersing branches contribute equally to the integral. Now we write − = δ n , with δ > 0, (so that δ = | | 1 n ). We then make the substitution x = δy. These then give: For > 0, the constant energy curves disperse rightward/leftward. The vertices are at (± 1 n , 0). The DOS integral in this case is: Where once again we have used the same set of substitutions to simplify the integrand. There are four four branches of the contour which contribute equally to the above integral. For = y 2 − x n with n even, the roles of positive and negative are interchanged in the above derivation.
For odd n the constant energy contours of = x n − y 2 disperse rightward for both 0 and > 0. When < 0, the vertex of the contour is (−| | 1 n , 0). The DOS is given by: Where once again we have scaled x to extract out the power law dependence. On the contrary, when > 0, the vertex is at (| | 1 n , 0) and the DOS is given by: The dispersion = y 2 − x n with n odd just has the role of positive and negative energies interchanged in the above. In fact = x n + y 2 is also equivalent to the above case. This is because, the coordinate transformation x → −x does not affect the DOS.
Thus, g( ) ∼ | | 1 n − 1 2 for the cuspoid catastrophes although the coefficient in front depends on the dispersion.

The umbilics and the rest
We first treat the hyperbolic umbilic and elliptic umbilic catastrophes and their higher order generalizations. These respectively take the form (x, y) = x 2 y ± y n for odd n. The energy contours of these are symmetric in the sense that the transformation y → −y takes us from the contour of to the contour of − . From this property, it is easy to see that the coefficients D + and D − are equal so that D + /D − = 1. Since we are ultimately interested only in their ratio, we will not bother evaluating them and just apply the appropriate scaling procedure to extract the exponent. To this end consider > 0 and set x = 1 2 − 1 2n and y = 1/n v. The Jacobian determinant for this transformation is: Using the scaling property of the delta function, viz δ(ax) = δ(x)/|a| the DOS integral becomes: The parabolic umbilic and its generalizations take the form x 2 y + y n for even n. We do not separately consider the case x 2 y − y n since it can be obtained from the former by the transformation y → −y followed by an overall sign change.
We proceed just as in the case of the cuspoids to compute D + /D − .

Appendix B: Symmetries and critical points
We now compute the effect of reflection (mirror) symmetry and rotational symmetry on the Jacobian. First consider reflection about the k x axis. This is represnted by: Now for a coordinate transformation φ the Jacobian at the origin transforms as Df (φ(0)) = Df (0)Dφ(0). Let k = (k x , k y ). Since f (k) = f (r k x · k), we have: (B2) This implies ∂f /∂k y = 0 at the origin. Thus, the direction perpendicular to the mirror has a vanishing component for the Jacobian. Now let ρ θ denote rotaion by θ in the plane. Applying the same procedure as above gives: The determinant of this homogeneous system is 4 sin 2 θ 2 , which is non zero for θ = π/3, π/2, π, 2π/3 which are generators of non-trivial lattice rotation symmetries. This implies the unique solution is ∂f /∂k x = ∂f /∂k y = 0, so that we have a critical point.

Appendix C: Quick review of catastrophe theory
The discussion in this section closely follows [20,21]. We first review a few definitions and conventions. If U is an open subset of R n , a function f : U → R is smooth if derivatives of all orders exist. From here on f will always denote a smooth function defined on U , an open subset of R n . We shall use map and coordinate transform interchangeably to refer to a bijective function φ : U → V , where U and V are open subsets of R n . The derivative of f at any point x 0 ∈ U , denoted by Df | x 0 : R n → R, is a linear transformation. When it acts on a vector, it gives the derivative of the function with respect to the vector; in particular, unit vectors give directional derivatives. Also, the function f (x 0 ) + Df | x 0 · (x − x 0 ) is the best linear approximation to f , near x 0 (i.e., in some neighborhood of x 0 ). These notions extend to higher derivatives as well. We note that if x is represented by a column vector, then Df | x 0 is represented by a row vector (that we refer to as the Jacobian), while the second derivative, denoted by D 2 f | x 0 is an n × n matrix called the Hessian.
x 0 is a critical point if Df | x 0 vanishes. Since, from now on we shall focus on critical points, we assume that the origin has been translated so that the critical point we are interested in occurs at the origin.

Diffeomorphism
exists an open neighborhood V of x 0 such that ψ restricted to V is a diffeomorphism. It can be shown that φ is a local diffeomorphism at x 0 if and only if det Dφ| x 0 = 0.

Equivalence of functions
Two funtions f and g are equivalent if there exists a diffeomorphism φ : U → V and a constant γ (called the shear term) such that g(x) = f (φ(x))+γ, ∀x ∈ U . We shall often refer to such a φ as a smoothly reversible coordinate transformation.

Morse and non-Morse critical points -corank
If the Hessian of f at origin is nondegenerate (i.e., it's rank equals n), then Morse lemma tells us that in some neighborhood of the origin, f is equivalent (in the above sense) to a quadratic form −y 2 1 − · · · − y 2 l + y 2 l+1 + · · · + y 2 n . A quadratic form having this structure is referred to as a Morse l-saddle. (Thus, a maximum is an n-saddle while a minimum is a 0-saddle). The critical point is itself referred to as a Morse critical point or an ordinary critical point. Another way to state Morse lemma is to say that near an ordinary critical point, there is a local coordinate system in which the function takes the form of a Morse l-saddle.
When the Hessian is degenerate and has rank r (and therefore corank n − r), then we can apply the splitting lemma which states that the function is equivalent near 0 to a function of form ±x 2 1 ± · · · ± x 2 r +f (x r+1 , ..., x n ). We refer to ±x 2 1 ± · · · ± x 2 r as the Morse part andf (x r+1 , ..., x n ) as the non-Morse part or residual singularity. The critical point is itself referred to as a non-Morse critical point or a higher order critical point. The corank of the function, denoted by cor[f ] is the number of zero eigenvalues of the Hessian.

Structural stability
We say f is structurally stable if for sufficiently small perturbations p, f is equivalent to f + p. We know that non-degeneracy of the Hessian is equivalent to the condition det D 2 f | 0 = 0. Since determinant is a continuous function, we expect that when we add a small perturbation p, det D 2 (f + p)| 0 is still non-zero in some neighborhood of 0, so that the critical point (which might have moved), is still Morse. In fact the converse is also true so that we have the following result: A critical point is structurally stable if and only if it is Morse.

Jets and jet spaces
The k-jet of a smooth function, denoted by j k f , is obtained by taking the Taylor series up to degree k (i.e., truncating terms of O(k + 1) and higher). The vector space J k n is the set of all polynomials of degree k in n variables, with zero constant term. (alternatively J k n = {j k f | all f : R n → R with f (0) = 0}).

Determinacy
We say that a two functions f and g are k-equivalent if j k f = j k g. A function f is k-determined if every function k-equivalent to it is also equivalent to it. In particular we have j k [j k f ] = j k f so that j k f ∼ f . This means that if f is k-determined, we can find a smoothly reversible coordinate transformation that maps f to j k f , effectively removing terms of O(k + 1) and higher from it's Taylor expansion.
The determinacy of a function, denoted by det[f ] is the smallest k for which f is k-determinate. If there is not finite such k then det[f ] = ∞. The determinacy of a function is preserved under a smoothly reversible coordinate transformation.

Codimension
We first define some polynomial spaces: E k n is the vector space of all polynomials in x 1 , · · · , x n of degree k. J k n is the subspace of E k n of polynomials with zero constant term. M k n is the subspace of all homogeneous polynomials of degree k.
The expression P k , where P is a polynomial, refers to the truncation of P to degree k.
The tangent space ∆ k (f ) to f is the subspace of J k n spanned by polynomials of the form Qj k ∂f It contains the directions in which j k f moves under smooth coordinate transformations. It can be shown that The codimension of f , denoted by cod(f ) is the codimension of ∆ k (f ) in J k n , for any k for which f is k-determinate. cod(f ) tells us the number of missing directions in the subspace in which j k f moves. Thus, these directions, which are really polynomial terms, can not be added or removed by smooth coordinate changes. The codimension of a function is preserved under a smoothly reversible coordinate transformation.

Winding number
We introduce the notion of winding number which counts the number of times the function changes sign along a closed contour around the origin. Since this is an integer, we expect that it does not change under smooth coordinate transformation.

Unfolding
An r-unfolding of f at 0 is a function F : We also refer to it as an r-parameter family through f . In this context, the x i are state variables while t i are control variables.
A d-unfolding F is induced from an r-unfolding F via three mappings defined in a region around the origin so that F (x, s) = F (y s (x), e(s)) + γ(s) where: e : R d → R r , (s 1 , · · · , s d ) → (e 1 (s), · · · , e r (s)), y : R n+d → R n , (x, s) → (y 1 (x, s), · · · , y n (x, s)), An r-unfolding of f at 0 is versal if all other unfoldings of f at 0 can be induced from it. If further r = cod(f ), then it is universal. Finally, we mention an important result: If f is k-determinate, then we can construct a universal unfolding for f by choosing a cobasis {p 1 (x), · · · , p r (x)} for ∆ k (f ) in J k n and defining: F (x 1 , · · · , x n , t 1 , · · · , t r ) = f (x)+t 1 p 1 (x)+· · ·+t r p r (x).

(C1)
Appendix D: Thom's theorem basic catastrophes The original Thom's theorem states that any r-parameter family of smooth functions R n → R for r 4 and any n 1 is typically structurally stable and equivalent to one of the following seven types: The fold has a universal unfolding: x 3 + t 1 x ± y 2 . The determinacy is 3 while codimension is 1, so that the fold can be thought of as two ordinary critical points merging into a higher order critical points. Only symmetry present in the standard form is reflection about y-axis (i.e., r y ).

b. Cusp
Any universal unfolding of the cusp is equivalent to: ±(x 4 + t 1 x 2 + t 2 x) ± y 2 . The determinacy is 4 and since the codimension is 2, under the action of the unfolding, three ordinary critical points merge to give a cusp. This could be preceeded by a situation wherein just two of the three critical points merge to give a fold and Morse critical point. This is prohibited if t 2 = 0, since in that case we have r x and ρ π symmetries in addition to r y .

c. Swallowtail
The swallowtail has a determinacy of 5 and is described by the unfolding: x 5 +t 1 x 3 +t 2 x 2 +t 3 x±y. With a codimension of 3 it corresponds to the merging of 4 ordinary critical points on a line. On the way to the swallowtail we can get cusps and folds. Similar to fold it has only r y as symmetry.

d. Butterfly
For the butterfly, the universal unfolding takes the form: ±(x 6 + t 1 x 4 + t 2 x 3 + t 3 t 2 + t 4 x) ± y 2 . This has the consequence that the 5 ordinary critical points lying on the x-axis merge to give higher order critical point(s). The unfolding can also exhibit the fold, cusp and swallowtail due to the merging of fewer than five critical points. When t 2 = 0 and t 4 = 0, it has r x and ρ π as symmetries in addition to r y .

e. Elliptic umbilic
The standard form of the unfolding for the elliptic umbilic is: x 3 −3xy 2 +t 1 (x 2 +y 2 )+t 2 x+t 3 y. With determinacy 3 and codimension 3 it corresponds to four critical points merging to give a monkey saddle. When t 2 = t 3 = 0, it possesses threefold rotation symmetry (ρ 2π/3 ) along with r y .

f. Hyperbolic umbilic
The hyperbolic umbilic is similar to the elliptic umbilic except for a sign difference; the unfolding is: x 3 +3xy 2 +t 1 x 2 + t 2 x + t 3 y. It too has determinacy 3 and codimension 3. However it does not possess the three-fold rotational symmetry of the monkey saddle: the only symmetry is r y when t 3 = 0.
We have illustrated Thom's theorem in two dimensions. For n dimensions, we simply add a Morse part ±y 2 1 ± · · · ± y 2 n−2 to the unfoldings given. Thoms theorem can be extended to r = 6 for arbitrary n and r = 7 for n = 2. This requires us to incorporate all of the catastrophes with codimension r.
This has been done in Tables II and III. In addition to these, we consider one other singularity which is allowed by π/4 rotation symmetry: h. X 9 This is actually a class of mutually inequivalent singularities of the form x 4 + 2cx 2 y 2 + y 4 . For |c| = 1 this is a higher order singularity with codimension 8 and determinacy 4. For two different parameters c and c , the singularities are equivalent if and only if c, c > −1 and c = (3 − c)/(1 + c). The function itself has a four-fold rotational symmetry and we will use a symmertry consistent unfolding: = x 4 + 2cx 2 y 2 + y 4 + t(x 2 + y 2 ).

Appendix E: Taylor expansions and catastrophes
A family which is smoothly equivalent to the standard universal unfolding of a catastrophe can have a Taylor expansion that does not resemble the unfolding to any order. In fact, it may not be obvious what catastrophe is equivalent to the given family by simply looking at the Taylor expansion truncated to any degree. This complication presents not just for families and catastrophes, but also for functions and standard forms of higher order singularities. Given a function with a higher order critical point at the origin, it is in general very hard to find the local diffeomorphism that transforms it to the standard form. So in practice, it becomes necessary to avoid finding explicit coordinate transformations in order to determine the type of higher order critical point.
As seen earlier, corank, codimension, determinacy, winding are preserved under a smoothly reversible coordinate transformation and the standard types are described uniquely by these numbers. Thus if we calculate these numbers for a given function, we can unamibguisly identify the type of higher order singularity it is equivalent to. It is fairly easy to compute the corank and winding of a function (see Appendix C 3 & C 8). The problem therefore reduces to computing the determinacy and codimension of the given function in an efficient way. This is achieved by the method given in Sec. VI A where we also apply the method to an example. We now explain the reason behind the method.

Explanation of the method
For any function f , the k th degree Taylor polynomial j k f moves in the space of all polynomials under a smooth coordinate transformation. For example, let j 3 f = k 3 x − k 2 y . Under k x → k x and k y → k y −k 2 x , we have j 3 f → k 3 x −k 2 y +2k 2 x k y . Assume that the critical point of the function occurs at the origin. By applying a family of origin preserving smooth coordinate transformations φ(t) parametrized by a continuous real number t, we can generate the orbits j k [f (φ(t))] of j k f in the polynomial space. (We assume that φ(0) is the identity coordinate transformation and refer to the original Taylor polynomial at the origin j k [f (φ(0))] as simply j k f ). The orbits are essentially curves γ(t) parametrized by t with γ(0) = j k f . We can define tangents to these curves at j k f : The space of such tangents for all possible origin preserving one parameter families is the tangent space at j k f . An important result in catastrophe theory guarantees the following: if the tangent space at j k f contains all homogeneous monomials of the form k j x k k−j y for j = 0, · · · , k, then terms of O(k + 1) can be removed from f by a smooth coordinate transformation. This ensures that f is k-determinate. Conversely, if f is (k − 1)-determinate, then the tangent space at j k f contains all the homogeneous monomials k j x k k−j y for j = 0, · · · , k. By determining the lowest k for which the tangent space at j k f contains k j x k k−j y for j = 0, · · · , k we can narrow down the determinacy to either k or k − 1. To this end, in the method in Sec. VI A we generate a set of polynomials {p i } that span the tangent space and try to take linear combinations of these polynomials to generate k j x k k−j y for j = 0, · · · , k.
If we also allow smooth coordinate transformations that move the origin, the tangent space gets enlarged. The codimension of f is then codimension of this enlarged tangent space at j k f in the space of all polynomials with zero constant term and degree k, for any k for which f is k determinate. Thus by finding a k for which f is k-determinate and generating a spanning set for the enlarged tangent space, we can find the codimension by finding the dimension of the tangent space and subtracting it from the dimension of the polynomial space. To find the dimension of the tangent space, we just list the spanning vectors as the rows of a matrix in the basis of monomials and compute the rank of this matrix.

Appendix F: Isometries of the plane
This section is mostly based on the treatment in [23]. Let P denote the two-dimensional plane. When an origin is chosen, P can be identified with R 2 so that any vector u can be written in terms of two real components as (u 1 , u 2 ) (here on, we will use P and R 2 interchangeably). For two vectors x = (x 1 , x 2 ) and y = (y 1 , y 2 ) let d(x, y) = (x 1 − x 2 ) 2 + (y 1 − y 2 ) 2 denote the euclidean distance between them. The euclidean distance can also be obtained from the dot product: x · y = x 1 y 1 + x 2 y 2 . An isometry of the plane is a function f : P → P which preserves the euclidean distance. i.e., d(f (x), f (y)) = d(x, y).
We denote a copy of R 2 as V, the group of all translations (it is a group under addition of vectors). It acts on P as follows: If a ∈ V, the translation by a, denoted by t a : P → P, simply adds a to the vectors in P: t a (x) = x + a. It is important to note that V has a fixed origin while for P, we can choose any point as the origin. The following are equivalent definitions of an orthogonal operator φ : R 2 → R 2 [23] (a) φ is an isometry that preserves the origin: φ(0) = 0, (b) φ preserves dot products: φ(x) · φ(y) = x · y, (c) The matrix of φ, say A, is such that it's transpose is it's inverse A T = A −1 .
Any orthogonal operator in two dimensions is either a rotation by some angle θ, which we denote by ρ θ or a reflection about some line through the origin. In the latter case, it can be uniquely decomposed as the composition of a rotation ρ θ and reflection about the x-axis (which we shall denote by r): φ = ρ θ r (see below). It is clear from above, that orthogonal operators act on P only after an origin has been chosen. Also, it is easily seen that translations and orthogonal operators are isometries. In fact, it can be shown that any isometry f can be uniquely decomposed as the composition of a translation and an orthogonal operator: f = t a φ [23]. Moreover, the set of all isometries of the plane is a group with the composition of functions as it's group law of composition. We denote this by M . It contains translations, rotations about points (not just the origin), reflections about lines and glide reflections.

Compositions of basic isometries
Using the matrix form of ρ θ and r, we can easily show that rρ θ = ρ −θ r. The operator ρ 2θ r is actually a reflection about the line through the origin which makes an angle θ with the x-axis. The easiest way to see this is to do a coordinate transformation U = ρ −θ which rotates the plane by −θ so that the line in consideration becomes the x-axis. Under this, an isometry φ transforms as φ → U φ U −1 . In particular we have ρ 2θ r → ρ −θ ρ 2θ rρ θ = ρ −θ ρ 2θ ρ −θ r = r. Thus, in the new coordinate system the isometry is the reflection about x-axis so that in the old coordinate system it is the reflection about the θ-line through the origin. Next we show that the composition of two reflections is a rotation: rρ θ rρ α = rrρ −θ ρ α = ρ α−θ . Now we investigate the composition of a translation and an orthogonal operator.
We first note the following: ρ θ t a = t ρ θ (a) ρ θ . To show that the composition of translation t a and rotation ρ θ is a rotation about some point b by θ, we need to show that the equation t a ρ θ = t b ρ θ t −b has a unique solution for b. (The right hand side amounts to shifting b to the origin, rotating the plane by θ and shifting back by b which is the same as rotating by θ about b). Since t b ρ θ t −b = t b−ρ θ (b) , the equation reduces to a = (1 − ρ θ ) b, which is a pair of linear equations in two unknowns (the two components of b). This has a unique solution if and only if det(1 − ρ θ ) = 0. But det(1 − ρ θ ) = 2 − 2 cos θ. For θ ∈ [0, 2π), only θ = 0 satisfies 2 − 2 cos θ = 0. Thus, for any θ = 0 (mod 2π), we have a unique non-zero b corresponding to each b = 0. This verifies the assertion that any composition of rotation and translations is a rotation by the same angle about some point.
Before we look at the composition of translations and reflections about lines through the origin (which take the form ρ θ r), we rotate the plane to make the reflection axis the x-axis. We first consider the composition of a y-translation with the x-reflection: t aĵ r = t a 2ĵ t a 2ĵ r = t a 2ĵ r t − a 2ĵ . But the latter expression really performs reflection about the line y = a/2 (since the operation shifts this line to the x-axis, performs x-reflection and then shifts backs the line). Thus, the composition of a reflection with translation by a vector perpendicular to the reflection axis amounts to simply shifting the reflection axis by half the translation vector. For a general translation vector, we split the vector into a parallel and perpendicular parts (with respect to the reflection axis). The perpendicular part shifts the reflection axis while the parallel part adds a glide term: t a r = t a t a ⊥ r = t a (t a ⊥ /2 r t −a ⊥ /2 ). Thus, a generic composition of a translation and a reflection is a glide reflection. As before, the composition of two glide reflections is a rotation.

Discrete subgroups of isometries
A subgroup G of M is discrete if it does not contain arbitrarily small rotations and translations. Mathematically, this means that there is a positive real number such that for any translation t a ∈ G, we have a < and for any rotation about some point, we have the angle of rotation θ < .
Given the unique decomposition of an isometry as f = t a φ, we can define a map π : M → O(2), π(t a φ) = φ. This is a group homomorphism. We now restrict π to G. The image of G under π, denoted byḠ is a discrete subgroup of O(2). This is called the point group. An element of point group need not be an element of G, so that in generalḠ is not a subgroup of G. (This can happen for instance if G contains only glide reflections and no pure reflections. But the point group will contain a pure reflection). Thus, the point group contains information about the orthogonal operations present in G, either as pure orthogonal operators or in combination with translations. We can show that the only discrete subgroups of O(2) are the finite subgroups: the cyclic subgroups C n and the dihedral groups D n . The cycilc subgroup C n of order n is generated by the rotation ρ θ with θ = 2π/n while the dihedral group D n , of order 2n is generated by ρ θ and r , where θ = 2π/n and r is a reflection about a line through the origin. Thus, for any discrete group of isometries G, the point groupḠ is C n or D n for a positive integer n.
There is one other group assosiated with a dicrete group of isometries: the lattice, denoted by L. It contains all the translations contained in G. In contrast to the point group, the lattice is indeed a subgroup of G. The lattice could be: (a) the zero group {0}, (b) the set of integer multiples of a non-zero vector a: L = {ma | m ∈ Z}, (c) the set of integer combinations of two independent vectors a and b: L = {ma + nb | m, n ∈ Z}.
There are precisely seventeen groups which have lattices of the form given in (c). They are known as the wallpaper groups and we shall restrict our interest to them. When the case (c) occurs, the point group is restricted to be C n or D n with n = 1, 2, 3, 4 or 6. This is known as crystallographic restriction. Another important property that ties up L andḠ is thatḠ preserves L. In other words, for any φ ∈Ḡ, φ(a) ∈ L ∀a ∈ L or φ(L) = L.