Qsymm: algorithmic symmetry finding and symmetric Hamiltonian generation

Symmetry is a guiding principle in physics that allows us to generalize conclusions between many physical systems. In the ongoing search for new topological phases of matter, symmetry plays a crucial role by protecting topological phases. We address two converse questions relevant to the symmetry classification of systems: is it possible to generate all possible single-body Hamiltonians compatible with a given symmetry group? Is it possible to find all the symmetries of a given family of Hamiltonians? We present numerically stable, deterministic polynomial time algorithms to solve both of these problems. Our treatment extends to all continuous or discrete symmetries of non-interacting lattice or continuum Hamiltonians. We implement the algorithms in the Qsymm Python package, and demonstrate their usefulness through applications in active research areas of condensed matter physics, including Majorana wires and Kekule graphene.


I. INTRODUCTION
A transformation that leaves a physical system invariant is called a symmetry, and such transformations have an ever-important role in modern physics.For example, symmetry breaking characterizes the classical theory of phase transitions, and the invariance of the speed of light between reference frames is a cornerstone of special relativity theory.In molecules and crystals, the symmetries of the constituent atomic orbitals determine the character of chemical bonds.The band theory of solids uses the translational invariance of a crystal structure to classify states into energy bands according to their crystal momentum, where the band structure is in turn constrained by the point group symmetries of the crystal.To describe such bands, model Hamiltonians based on tightbinding approximations [1][2][3][4][5] or k•p perturbation theory 6,7 are typically constructed by fitting a generic Hamiltonian allowed by symmetry to match experimental data or first principles-calculations.
Recent studies focused on the role of symmetry in protecting topological phases. 8,94][15][16][17][18][19] Some of these phases are stable against disorder that preserves the symmetry only on average, 20 leading to a richer structure of symmetry-protected topological phases.Analysis of newly proposed symmetryprotected topological phases is often done using minimal models, such as tight-binding Hamiltonians with shortrange hoppings, or continuum Hamiltonians valid near high symmetry points in the Brillouin zone.Although these models are usually easy to solve, they are prone to having higher symmetry than intended.With the plethora of available symmetry groups, it is a nontrivial task to construct models that possess the stated, but only the stated symmetries, or to decide the complete symmetry group of a given Hamiltonian.In this paper we present an algorithm to generate all Hamiltonians that respect given symmetries, using an approach similar to Ref. 21.In addition, we develop a dual algorithm to find all symmetries of a family of Hamiltonians (Fig. 1).Our framework is applicable to all non-interacting, finite or translation invariant lattice or k • p Hamiltonians.We treat all possible symmetries, including continuous unitary symmetry groups, continuous spatial rotations, space groups, discrete onsite symmetries (such as time reversal, particle-hole and chiral symmetries), and arbitrary combinations of these.arXiv:1806.08363v1[cond-mat.str-el]21 Jun 2018 Besides static fermionic systems, our methods are also applicable to band structures in photonic crystals, 22,23 magnon spectra, classical mechanical 24,25 and electronic systems, 26 and driven Floquet systems. 27he paper is structured as follows: first we review the general symmetry structure of single-particle Hamiltonians.Then we present our algorithm to find Hamiltonians with such symmetries.After that we review the symmetry finding algorithm, which, by factoring out onsite unitary symmetries, makes finding all other symmetries more efficient and guaranteed.We implement our algorithms in the Qsymm Python package. 28,29Finally we provide a set of examples where our method was used on problems in active areas of research.We show that Majorana wire devices may be protected against bandtilting by a magnetic symmetry, and double Dirac cones in Kekule distorted graphene are protected by point group and sublattice symmetry.We also construct model Hamiltonians for transition metal dichalcogenides and distorted spin-orbit coupled SnTe.

II. HAMILTONIAN FAMILIES AND SYMMETRIES A. Continuum and tight-binding Hamiltonians
We focus on non-interacting Hamiltonians.The quadratic Hamilton operator of a finite (zerodimensional) system can be written as where H is a Hermitian matrix and âi are creation or annihilation operators.We do not make any assumptions about bosonic or fermionic nature of these operators and also allow terms with two creation or two annihilation operators, facilitating the study of superconducting Bogoliubov-de Gennes Hamiltonians.In the framework we use, all the details of the system are encoded in constraints on the matrix part H, which is the focus of our study.
Besides finite systems, we also address systems with ddimensional translation invariance.The associated conserved quantity is the (lattice) momentum k, which takes values in R d for continuous translations, and in the ddimensional Brillouin zone for discrete translations.Effective continuum Hamiltonians (k • p models) also arise as the long-wavelength limit of lattice Hamiltonians.The conservation of k allows to decompose the single-particle Hilbert space into independent sectors corresponding to each k, such that Ĥ does not mix sectors, i.e.
In the rest of this work we focus on analyzing the matrixvalued Hamiltonian H(k) with matrix elements H ij (k).
A tight-binding Hamiltonian acts on a Hilbert space consisting of basis orbitals in a single translational unit cell, and has the general form The hopping vectors δ connect sites on the lattice, with the matrices of hopping amplitudes h δ .The k•p Hamiltonian provides an accurate continuum approximation near a point in the Brillouin zone, typically a high symmetry point.It is a polynomial in momenta k i , and has the form where is a monomial in the multi-index notation with n = (n 1 , n 2 , . ..), and h † n = h n .Typical methods to construct k • p Hamiltonians start with a series expansion of a more complete lattice model from e.g.[32]

B. Hamiltonian families
A set of symmetries only defines a Hamiltonian family, as opposed to one single Hamiltonian.A Hamiltonian family is the linear space of Hamiltonians with arbitrary real coefficients c n , and basis vectors Here h αn are constant matrices of identical size, and f α (k) are linearly independent scalar functions.In the rest of this work, whenever referring to a Hamiltonian, we mean a Hamiltonian family.A Hamiltonian family is also the only useful starting point to analyzing the symmetry content of Hamiltonians.For a zero-dimensional Hamiltonian, the symmetry group is always given by independent unitary transformations in each degenerate eigensubspace.This group, however, provides no insight beyond the degeneracies of the levels.In a family of continuum k • p Hamiltonians, f α (k) is a monomial, and in a tight-binding Hamiltonian, a phase factor e ik•δ with δ a hopping vector.
There is a natural inner product in the space of Hamiltonians (5).We define the product of On the matrix part, the Frobenius inner product is given by A, B F = Tr A † B .For the inner product of the k-dependent prefactors, we use the Bombieri inner product 33 for polynomials, such that different monomials are orthogonal, and e ik•a , e ik•b = δ a,b for phase factors.Both of these inner products on the function spaces are invariant under the isometries of kspace, and therefore all symmetry actions we consider in this work are (anti)unitary with respect to this inner product.This structure of the space of Hamiltonians also justifies our use of single exponentials and single monomials as the expansion basis.

C. Symmetry constraints on Hamiltonian families
We adopt the active view of symmetry action g on the Hamiltonian: g(H) represents a transformed Hamiltonian, such that the matrix elements between rotated wave functions g (|ψ(k) ) are identical.In other words, a Hamiltonian has a symmetry if g leaves the Hamiltonian invariant, A general unitary symmetry g acts on a Hamiltonian H(k) as and a general antiunitary symmetry as Here the orthogonal matrix R g is the real space action, and the unitary matrix U g is the Hilbert space action of g.We include the overall ± sign to treat antisymmetries-symmetries that reverse the sign of energy-on an equal footing.We restrict our considerations to a constant U g , however, in the real space basis (see Section III B), any space group operator may only contain an overall k-dependent phase factor, which cancels in the previous equations.Substituting Eq. ( 7) into ( 6), we rewrite it in a form linear in the symmetry action: Here the symmetry action S is S = U if unitary and S = U K if antiunitary, with K the real space complex conjugation operator: We apply the symmetry constraint (8) to the Hamiltonian family (5a).The spatial action of a symmetry is a rotation in the space of f α , such that f α (±Rk) = β γ β α f β (k) with γ β α known for given R and f α .Substituting this yields Since f α are linearly independent functions, the matrix coefficients in the parentheses must vanish for every α, resulting in the system of equations When symmetries form continuous Lie groups, it is advantageous to use the symmetry generators instead of the group elements.Consider a one parameter family of transformations g φ which, for a fixed φ act as a unitary symmetry in the above, and let g 0 be the identity.We define the action of the generator g through Substituting (7) and using that U g φ = e −iφL with L = L † and R g φ = e −iφM with M = −M T = M † , we find g φ (H) = H for every φ is equivalent to g (H) = 0. Tightbinding Hamiltonians cannot be invariant under continuous rotations of space, such that M = 0 and the symmetry constraint simplifies to where L is a local conserved quantity.Finally, if the Hamiltonian is a polynomial in k, continuous rotation invariance is also possible.With

A. Constraining Hamiltonian families
Given a symmetry S and a Hamiltonian family (5), we wish to find the subfamily of Hamiltonians that is invariant under the symmetry transformation (8).The symmetry constraint on the Hamiltonian family is a system of homogeneous linear equations for the coefficients c n [see Eqs.(10), (13), and (14)].We find the space of solutions numerically using singular value decomposition or sparse eigendecomposition, which gives the subfamily of the original Hamiltonian family (5) that satisfies the symmetry.Imposing additional symmetry constraints on the family yields further linear equations that are identical to Eq. (10) in form.We provide an implementation of this algorithm in the Qsymm Python package.
The constraining algorithm allows to generate all possible tight-binding or k•p Hamiltonians that satisfy symmetry constraints, by applying the algorithm to the most general representative Hamiltonian family for the system at hand.As an illustration, we reproduce the family of two-band k•p Hamiltonians of Ref. 30 for the surface dispersion of the topological insulator Bi 2 Te 3 .Our starting point is the family of all 2 × 2 k • p Hamiltonians up to third order in the momentum k = (k x , k y ).Expanding the matrix part in terms of the identity and Pauli matrices σ 0,x,y,z , the general family consists of 40 basis vectors and is given by with n = (j, α x , α y ) = (j, α).To obtain the surface dispersion Hamiltonian, we constrain (15) with timereversal symmetry (T = iσ y K), and the point group symmetries of the crystal, namely three-fold rotation (C 3 = e −iπσz/3 ), and mirror symmetry in x (M = iσ x ). 30ubstituting the three symmetries and the family ( 15) into (10) yields a homogeneous system of 120 linear equations for the 40 coefficients c n .The null space of the linear system is the subfamily of Hamiltonians that satisfy the symmetry constraints: with c n ∈ R, which matches the Hamiltonian of Ref. 30.
Here, k 2 = k 2 x +k 2 y , and we have relabelled the coefficients c n for clarity.

B. Generating lattice Hamiltonians by symmetrization
Lattice models often contain multiple sites per unit cell, but only a small number of bonds.In this case the previous approach of generating all possible terms and constraining them is inefficient due to the large dimension of the null-space.On the other hand, all the hopping terms on symmetry equivalent bonds are completely determined by the hopping on one of these bonds.This allows us introduce a symmetrization strategy to generate all symmetry-constrained lattice Hamiltonians with hoppings of limited range.
To treat arbitrary space group symmetries of general crystal structures, we consider tight-binding Hamiltonians in the real space basis that preserves information on the coordinates of the basis orbitals. 34,35Up to a normalization factor, the Bloch basis functions are given by χ al k = R e ik(R+ra) φ al R , where a indexes the sites in the unit cell, l ∈ [1, . . ., n a orbs.] indexes the orbitals on the site, r a is the real space position of the site and R runs over all lattice vectors.In this basis, the hopping terms in the Hamiltonian acquire a phase factor corresponding to the true real space separation of the sites they connect, as opposed to the separation of the unit cells to which the sites belong.A hopping between site a at r a and site b at r b = r a + δ ab enters as a term e ikδ ab h ab δ ab + h.c.where we suppressed the orbital indices of the matrix h ab δ ab .Onsite terms have δ aa = 0.The main advantage of using this gauge is that the form of the Hamiltonian is independent of the choice of the real space origin and the shape of the unit cell.As a consequence, nonsymmorphic symmetry operations only acquire k-dependence in the form of an overall phase factor. 35In the simplest case of a single site per unit cell, δ ab are lattice vectors.
We start from a small set of terms for every symmetry unique bond δ ab of the form with h ab nδ ab spanning all n a orbs.× n b orbs.matrices that are invariant under the continuous onsite symmetry group.We symmetrize these with respect to the discrete point group G, i.e.
where g(H) is the symmetry transformed image of H under the transformation g [see Eq. ( 6)] and |G| is the number of elements in G.Because the sum over g ∈ G can be replaced by a sum over (hg) ∈ G, h(H s ) = H s for all h ∈ G.In addition, H s is exactly the projection of H onto the space of symmetric Hamiltonians.The symmetrized terms span all symmetry allowed Hamiltonians with the prescribed hopping vectors.This space of Hamiltonians is generally overcomplete, we find a minimal set of terms spanning the space using standard linear algebra techniques.
As an example, consider graphene with one spinless orbital per site.A three-fold rotation around a site maps both sublattices onto themselves, so the unitary part of the symmetry action is 1 2×2 .Let δ = (a 0 , 0) be a vector connecting nearest neighbors, and the corresponding hopping term which is Hermitian and only connects the two sublattices.
After symmetrization with respect to the full hexagonal group we obtain the well known minimal tight-binding model for graphene: 2 ky .

IV. SYMMETRY FINDING
Unlike finding a family of symmetric Hamiltonians, that amounts to solving a linear system, finding the sym-metries of a Hamiltonian family is more involved.We first focus on finite (zero-dimensional) systems and show that the unitary symmetry group generally admits a continuous Lie group structure.Next we present an algorithm to find the unitary symmetries, and rewrite the Hamiltonian family in the symmetry-adapted basis.In this basis the Hamiltonian takes a block diagonal form, where the blocks are guaranteed to have no unitary symmetries, hence we call these blocks reduced Hamiltonians.Factoring out the unitary symmetries this way simplifies finding the discrete (anti)unitary (anti)symmetries, see Fig. 2.After generalizing these methods to onsite symmetries of translation invariant systems in arbitrary dimensions, we finally include real space rotation symmetries.

A. Structure of the onsite unitary symmetry group
Assume a unitary symmetry operator U commutes with a family of finite Hamiltonians: Any unitary is expressible as the exponential of a Hermitian operator which is unique if we restrict the spectrum σ(L) ⊂ [0, 2π) (this condition is equivalent to choosing a branch cut for the logarithm in L = i log U ). Defining U and L this way ensures that they have exactly the same eigensubspaces.Since U and H n commute, they also share eigensubspaces, and Eq. ( 21) is equivalent to Therefore the Hamiltonian family has a continuous family of unitary symmetries U (φ) = e −iφL that all commute with H n .Because a single unitary symmetry defines a continuous family that is connected to the identity, the full group of unitary symmetries must form a single connected component G 0 .This is uniquely specified by the space of conserved quantities L ∈ g, the Lie algebra of the Lie group G 0 .
Consider for example a system consisting of a number of spinful orbitals, which is invariant under the set of unitary spin-flip operators U i = iσ i ⊗ 1, where σ i are the Pauli matrices and the identity acts on the space of orbitals.Taking the logarithms of these operators, we find that they are associated with the Hermitian conserved spins L i = σ i ⊗ 1.The Lie group generated by these conserved quantities is SU (2) acting in spin space.Therefore the generic Hamiltonian of such a system assumes the form H = 1 2×2 ⊗ H r , where the reduced Hamiltonian H r acts only on the space of orbitals and the identity acts on spin space.In the basis where spin up and spin down states are grouped together, the original Hamiltonian takes the block diagonal form with two identical blocks H = H r ⊕ H r , reducing the problem to spinless fermions.
Consider the same system, but with higher spins on every site instead.Let the conserved spins be L i = J i ⊗ 1 where We therefore find that the symmetry group is in fact larger than the one we started with, forming a full unitary group.The above result may sound surprising on physical grounds, considering that many well-studied models (e.g. the transverse field Ising model) only have discrete onsite symmetry groups.We emphasize that this result (and much of what follows) is specific to single-particle systems, and does not directly apply to onsite symmetries on the many-particle Fock space.In the single-particle case the full Hilbert space is the direct sum of the local Hilbert spaces and therefore an onsite symmetry is the direct sum of local unitaries.The many-particle Hilbert space is a direct product, and onsite unitary symmetries take a direct product form.The associated L is generally not a sum of local terms, and does not correspond to a local conserved quantity.The above argument also fails when considering spatial symmetries, because in general the logarithm of a locality-preserving operator (an operator that maps a state with localized support to one with localized support) mixes degrees of freedom that are far apart.
In the rest of this subsection we prove, using the theory of Lie groups, that the unitary symmetry group G 0 is a direct product of unitary groups U (N ) in any finite system.We then show the existence of the symmetryadapted basis, where both the conserved quantities and the Hamiltonian take a simple form, and derive properties of reduced Hamiltonians.The reader not interested in mathematical proofs may skip to the next subsection where the algorithm for finding unitary symmetries is discussed.
The unitary symmetry group G 0 is a subgroup of the full unitary group on the Hilbert space H, G 0 ≤ U (dim H).G 0 is a connected and compact matrix Lie group, which means that all of its finite-dimensional representations are completely reducible. 36,37The Lie group G 0 is generated by all the generators in its Lie algebra L ∈ g for which [L, H] = 0, the Lie algebra g is also completely reducible.
Reducing the representation amounts to splitting the Hilbert space H into a direct sum of irreducible subspaces Each of these subspaces is invariant under the symmetry action and contains no invariant subspace.V (i) j transforms according to irreducible representation (irrep) i, and irrep i has multiplicity n i .We denote the union of all irreducible subspaces belonging to irrep i as In a symmetry-adapted basis, every symmetry generator takes the same block diagonal form of irreducible representations: (25)   where L (i) is the representation of L in the i-th irrep, each acting in a corresponding irreducible subspace V (i) j of dimension d i .Irreps j and k are equivalent if there exists a unitary transformation W such that L (j) = W L (k) W −1 ∀L ∈ g.This guarantees that there is a basis where all operators have exactly the same representation in every equivalent irreducible subspace.
In this basis, the Hamiltonian also takes a simple block form.By Schur's Lemma, blocks of H between irreducible subspaces that transform according to different irreps are zero, and blocks between irreducible subspaces with identical irreps are proportional to the identity: where the reduced Hamiltonians H i are n i ×n i Hermitian matrices.The reduced Hamiltonians H i cannot have any nontrivial unitary symmetries.To prove that, assume that H 1 has a conserved quantity L such that [H 1 , L] = 0.It implies that (1 d1×d1 ⊗ L) ⊕ 0 ⊕ . . .commutes with the full Hamiltonian, which is incompatible with the unique decomposition to irreducible subspaces, except the trivial case of L ∝ 1.
It is apparent from (26) that the symmetry group of H is a product of full unitary groups acting independently on each block where the symmetry generators have the form (25), with L (i) ∈ u(d i ) independently running over all d i × d i Hermitian matrices, and u(d i ) the Lie algebra of U (d i ).Because the reduction to irreducible subspaces is unique, this is the full group of unitary symmetries.The center of the group Z(G 0 ) is formed by the abelian U (1) subgroups generated by the set of projectors on each block, i.e. generators where one of the L (i) is the identity and the others vanish.
To have a basis-independent characterization of the center of the Lie algebra, we compute the structure constants f αβ γ defined by Using these we define the Killing form It can be shown 37 that the null-space of the Killing form is exactly the center of the Lie algebra, i.e. if a vector l is a solution of β K αβ l β = 0 then α l α L α commutes with every operator in g.

B. Finding the unitary symmetry group
We are now ready to define the algorithm of finding the unitary symmetry group and constructing the reduced Hamiltonians for a given family of Hamiltonians.First we find all symmetry generators L α as the linearly independent solutions of This is a system of linear equations for the unknown components of L α , which we solve using the same methods we used for constraining Hamiltonians.After computing the Killing form K we find all linearly independent solutions of for l.Operators of the form α l α L α are the basis of conserved quantities that commute with every other conserved quantity.We simultaneously diagonalize all of these (see Appendix A) to find the simultaneous eigensubspaces V (i) (24).
We then find the generators L (i) of the SU (d i ) symmetry group of each block.To do so, we project the Hamiltonian onto V (i) using the projector P i , which is an orthonormal set of column vectors, and solve [L (i) , P † i HP i ] = 0 and and Tr L (i) = 0, (32) to find the d 2 i − 1 linearly independent solutions for L (i) .The final step is finding a basis within V (i) that gives the tensor product structure of ( 25) and ( 26) (see Appendix B).We use this basis and the resulting reduced Hamiltonians in the following.
These symmetries also form continuous families, because combining them with any onsite unitary symmetry also results in a discrete onsite symmetry of the same type.Because there is no continuous way to interpolate between unitary and antiunitary symmetries or between symmetries and antisymmetries, each type forms a disconnected component of the onsite symmetry group (see Fig. 2).To find one representative of each type of the discrete onsite symmetries, we utilize the symmetryadapted basis and reduced the Hamiltonian found in the previous subsection.The reduced Hamiltonians have no residual symmetries, which makes the discrete onsite symmetries unique and allows us to efficiently find them.
We start with time reversal symmetries of finite (zerodimensional) systems.T is a time reversal if Writing T = U K with unitary U and complex conjugation K, we obtain This is a nonlinear system of equations, and it is in general hard to solve.We show that using the reduced Hamiltonian simplifies it to a linear problem.We first consider a Hamiltonian that has one set of identical irreducible subspaces, i.e.H = 1 ⊗ H 1 .By (25), all conserved quantities have the form L = L (1) ⊗ 1, and span the full space of Hermitian matrices on the first Hilbert space of the tensor product.If L is a conserved quantity, so is T LT −1 , which implies U gU † = g.Therefore the unitary part of T is a direct product of two unitaries, U = V ⊗ W (see Appendix C), with V an arbitrary unitary matrix.Because V = 1 commutes with all unitary symmetries, we call T = 1 ⊗ W K the canonical time reversal symmetry.Due to the tensor product structure of T , (34a) reduces to Importantly H 1 has no unitary symmetries.Any nonzero solution W of ( 35) has ker W = 0: otherwise either ker W is an invariant subspace of H * 1 or ker H * 1 is nonzero, both incompatible with H 1 having no unitary symmetries.Considering two solutions W and W of (35) we find Because H 1 has no unitary symmetries W W † ∝ 1, which proves that any solution of ( 35) is unique and unitary up to a constant factor.In other words, any normalized solution of (35) automatically satisfies (34b).
As an example consider the reduced family of Hamiltonians This family has no residual symmetry, and solving (35)  we find that W has to commute with both σ x and σ z , so W ∝ σ 0 .Therefore H is invariant under T = σ 0 K, which is unique up to a phase factor.As a second example consider the reduced family with τ i the Pauli matrices.The condition (35) implies that W commutes with τ x σ 0 , τ z σ 0 , τ y σ y and anticommutes with τ y σ x , τ y σ z .The only solution is W ∝ τ 0 σ y , and T = τ 0 σ y K up to a phase factor.In the general case of multiple irreps, we find (see Appendix D for details) that a time reversal can only mix subspaces V (i) and V (j) if they correspond to irreps of the same dimensionality and multiplicity.The block structure of U has to be symmetric, it can only exchange subspaces pairwise or leave subspaces invariant.In order to find a time reversal, we iterate over all symmetric permutations of compatible subspaces, and check if a time reversal exists with the given block structure.Specifically we consider two reduced blocks of the Hamiltonian that are interchanged Because H i and H j have no unitary symmetry, the square of time reversal must have the form Therefore, following Appendix D we search for a time reversal of the form The relation T r H r = H r T r then reduces to This is the key result of this section, and it is a generalization of (35).Following the same reasoning as before, we conclude that any nonzero solution W ij of ( 42) is unique and unitary up to a constant factor.
As an example of this case consider the Hamiltonian family Solving (42) for W 12 we find that [W 12 , τ i ] = 0 for i = x, y, z, so W 12 ∝ τ 0 and T = τ 0 σ x K up to a phase factor.We also confirm that (42) does not have any solutions for W 11 or W 22 , so the block form of T is unique.Likewise, the canonical unitary and antiunitary antisymmetries act as a tensor product U ij = 1⊗W ij in each block, and either leave subspaces invariant or pairwise exchange compatible subspaces.The results analogous to Eq. ( 42) for unitary and antiunitary antisymmetries are respectively:

D. Onsite symmetries of k-dependent Hamiltonians
The above methods extend to the onsite symmetries of k-space Hamiltonians of arbitrary dimensions.An onsite unitary symmetry acts locally in k-space and is independent of k.Given a family of Hamiltonians H(k, α), we treat linearly independent functions of k as additional free parameters and apply the methods of section IV B.
We now turn to time reversal symmetry, which requires special treatment because it transforms k to −k.Because H(−k) is a reparametrization of the same Hamiltonian family, it is reduced if H(k) is reduced.The generalization of (42) to the k-dependent case is By the same argument as before, separating the Hamiltonian to irreducible blocks guarantees that the nonzero solutions are unique and unitary up to constant factors.The analogous results are true for particle-hole and chiral symmetry.

E. Point group symmetries
The point group of a crystal is always a subgroup of the finite point group of its Bravais lattice.Therefore we search for point group symmetries by enumerating possible real space rotations R g , and applying methods similar to the previous subsections to find whether it is a symmetry with appropriate U g .
Like discrete onsite symmetries, point group symmetries may be combined with onsite unitaries, forming continuous families.This ambiguity is again removed by using the reduced Hamiltonian.The analogous result to (42) for point group symmetries is (U g ) ij = 1⊗W ij where the blocks W ij satisfy Here W only has one nonzero block per row and column, and nonzero blocks only between compatible subspaces.If the order of the symmetry is greater than 2, permutations which are not symmetric are also possible.Because both H j (k) and H i (R g k) are reduced, the nonzero solution for W ij is unique and unitary up to normalization and a phase factor.With the knowledge of the full point group, the arbitrary phase factors appearing in W ij may be fixed such that the U g form a (double)group representation of the point group.
A similar argument applies to the case of antiunitary point group symmetries (magnetic group symmetries) and antisymmetries that involve spatial transformations.The analogous equations for unitary antisymmetries antiunitary (anti)symmetries are respectively:

F. Continuous rotations
To find continuous rotation symmetries of k • p Hamiltonians we utilize the symmetry-adapted basis of the onsite unitary symmetries again.Unlike discrete symmetries, the unitary action of a continuous symmetry cannot mix different blocks, because it continuously deforms to the identity.Therefore we treat reduced Hamiltonians H i separately.
In order to find a continuous symmetry generator g as defined in sec.II C, we simultaneously solve for every j, with constraints L (j) = (L (j) ) † , Tr L (j) = 0 and M = −M T = M † .We then expand H j (k) in a basis of monomials, and reduce (48) to a system of linear equations for the entries of L (j) and M .

V. APPLICATIONS
We implemented the symmetric Hamiltonian generator and symmetry finder algorithms of the previous sections in the Qsymm Python package. 28,29We provide an interface to define symbolic expressions of symmetries and Hamiltonian families using Sympy 38 and Kwant. 39Efficient solving of large systems of linear equations is achieved using ARPACK 40 (bundled for Python by Scipy 41 ).We provide the the source code with instructive examples as a software repository. 28,29The following examples illustrate how the algorithms were used to solve open research problems in condensed matter physics.We also provide the Jupyter notebooks 42 generating these results.

A. Symmetries of Majorana wire
An early version of our symmetry finding algorithm was used in Ref. [43] to find the symmetries of a superconducting nanowire in an external magnetic field.The analysis revealed unexpected symmetries of the model system in certain geometries that prevent band tilting and closing of the topological gap.Here we revisit these results.
The system under consideration is an infinite nanowire along the x axis with a semiconducting core and a superconducting shell covering some of the surface.In the presence of a magnetic field along the wire and a normal electric field, the wire undergoes a topological phase transition.This is marked by the gap closing and reopening, with Majorana zero modes appearing at each end of a finite wire segment. 44A component of the external magnetic field normal to the wire axis breaks the symmetry of the band structure (E(k) = E(−k)), leading to tilting of the bands and closing of the superconducting gap at finite momentum.
Reference [43] studied a tight-binding Hamiltonian of the wire and observed that depending on the geometry the band tilting may or may not occur.Applying the symmetry finder algorithm to wires with small crosssections (Fig. 3), we identify the key difference in symmetry.If the wire geometry has a mirror plane including its axis (M y ) as in Fig. 3(a), with external fields E z, and B x, we find that the symmetry group consists of 8 elements.The three generators of this group are particle-hole symmetry (P), a mirror plane perpendicular to the wire axis (M x ) and the combination of M y with time reversal T .This last symmetry M y T includes both a spatial transformation and time reversal, and is easily overlooked.This operator can be further combined with FIG. 3. Sketch of two possible geometries for the nanowire (yellow) with a superconducting shell (red).In (a), the geometry respects the mirror symmetry My, which gives rise to a chiral symmetry.In (b) however, the superconducting shell breaks the mirror symmetry, and hence the chiral symmetry is also absent.Although the sketches show finite segments, the wire is translationally invariant along its axis x in both cases.
particle-hole to result in a chiral symmetry C = M y T P, as pointed out in the earlier work.A nanowire enhanced by such an effective time reversal symmetry belongs to class BDI and supports multiple Majorana modes at its end. 45ymmetries M x and M y T require E(k x ) = E(−k x ) and prevent band tilting.Adding nonzero B z to the magnetic field breaks M x , but preserves M y T , still forbidding band tilting.Further reducing the symmetry by moving the position of the superconducting cover as in Fig. 3(b), or by applying B y breaks all symmetries relating k x to −k x and enables band tilting.

B. Kekule distortion in graphene
The Kekule distortion of graphene is a periodic pattern of weak and strong bonds tripling the size of the unit cell.In the Kekule-O pattern, weak bonds around plaquettes resemble benzene rings, and in Kekule-Y, strong bonds form Y shapes around sites (Fig. 4).After folding back the Brillouin zone, the K and K points are both mapped to the Γ point.A suitable mass term can now open a gap in the band structure.A recent work 46 reported that unlike the Kekule-O distortion, 47 the Kekule-Y distortion does not open a gap: instead it preserves a double Dirac cone at the Brillouin zone center.Using our algorithms we identify the symmetries protecting this double Dirac cone.First we find all the symmetries of the effective 4-band k • p model of Kekule-Y: with v 1 , v 2 band structure parameters.We find that it is symmetric under the full hexagonal point group D 6 (in fact the linearized model has a continuous rotation symmetry), time reversal and sublattice symmetry, which results from the bipartite nature of the honeycomb lattice.Next, we systematically generate all subgroups of this symmetry group and the corresponding symmetryallowed 4-band k • p Hamiltonians.We find that at least one antisymmetry is required to forbid a constant mass term that would open a gap at Γ.A minimal subgroup protecting the double Dirac cone is generated by sublattice symmetry and three-fold rotations.Removing sublattice symmetry, even while keeping the full D 6 point group, removes the protection of the double Dirac cone.Sublattice symmetry is broken by adding second neighbour hopping, or a staggered onsite potential compatible with the lattice symmetries in the tight-binding model.Symmetry finding shows that the effective model of Kekule-O, with v, ∆ ∈ R, has the same symmetry group structure.However, the mass term σ z τ x is allowed even in the presence of the full symmetry group.The key difference is the unitary action of rotations: the generator of continuous rotations is σ z + τ z in the Kekule-Y case, while it is σ z in the Kekule-O case.Sublattice symmetry is C = σ z τ z in both cases.Therefore no constant matrix can simultaneously anticommute with C and commute with the rotation generator in Kekule-Y, while a mass term is allowed in Kekule-O.This difference in the transformation properties stems from the different Wyckoff positions of the lattice sites.In Kekule-Y the three-fold rotation centers lie on lattice sites, while in Kekule-O the three-fold rotation centers lie at centers of hexagonal plaquettes.Using the tightbinding Hamiltonian generator, we confirm that the representation of 3-fold rotations in the low energy subspace at the center of the Brillouin zone is different for the Kekule-O and Kekule-Y systems.

C. k • p model of distorted SnTe
The cubic rocksalt material SnTe is the first example of topological crystalline insulators. 48Recently, using our method, Ref. 49 proposed that structural distortions can give rise to Weyl and nodal-line semimetal phases in the same material.Here we review these results.
The band gap of the cubic phase is smallest at the L point in the face-centered cubic Brillouin zone.We construct an effective k • p model expanded up to second order in k around L. The model has two orbital degrees of freedom, spanned by p orbitals on Sn and Te sites.The initial symmetry group of the L point is D 3d which is generated by inversion I, a three-fold rotation C 3 about the ΓL axis, and a reflection M x about the mirror plane containing Γ and two L points.Furthermore, the model should be invariant under time reversal T .The corresponding representations of the symmetry operators, listing the (anti)unitary action first and the k-space action second, are as follows, where φ = 2π 3 , and Pauli matrices σ i and s i act on orbital and spin degrees of freedom respectively.k is the momentum vector measured from an L point in a coordinate system where the z axis is aligned with ΓL (e.g.[111]) and the x axis is normal to a mirror plane (e.g. [110]).The σ z in the unitary action of inversion is a result of considering an L point, because the inversion center is one of the sites in the unit cell of the rocksalt structure, the other site is translated by a lattice vector under inversion and acquires a phase factor at nonzero momentum.
Applying the k•p Hamiltonian generator algorithm we find 8 symmetry-allowed terms.Ignoring the 3 terms that are proportional to the identity and do not influence band topology, we obtain the following Hamiltonian family: Breaking the three-fold rotational symmetry results in 8 new terms, 6 of which are not proportional to the identity: Further breaking inversion symmetry produces 22 additional terms, none of which is proportional to the identity.

D. Three-orbital tight-binding model for monolayer transition metal dichalcogenides
Monolayers of transition metal dichalcogenides MX 2 (M = Mo, W; X = S, Se, Te) are promising materials for use in electronics and optoelectronics. 50When doped, the MX x monolayers also become superconducting. 51In the 1H stacking, a monolayer consists of a layer of transition metal atoms M sandwiched between two layers of chalcogen atoms X.Each layer separately is a triangular Bravais lattice, with the X atoms in the top and bottom layers projecting onto the same position in the plane of M atoms, forming an overall honeycomb lattice.In the normal state, the monolayer is a semiconductor, with conduction and valence band edges at the corners of the hexagonal Brillouin zone ±K.Using that the wave functions at the band edges are predominantly composed of d-orbitals on the M atoms, Liu et al. [4] developed a three-orbital tight-binding model with nearest neighbor hopping.This model satisfies the symmetry group of the monolayer, and has band edges near ±K.Here, we reproduce their spinless tight-binding model using our algorithm for symmetric Hamiltonian generation.
The tight-binding basis consists of three d-orbitals on the M atom, namely Because the model does not include any orbitals on the X atoms, it has a triangular lattice, with lattice vectors a 1 = x and a 2 = (x + √ 3ŷ)/2.The symmetry generators are time reversal symmetry T , mirror symmetry M x , and three-fold rotation in the monolayer plane C 3 which are represented in the tight-binding basis as with φ = 2π 3 .Employing the symmetrization strategy for lattice Hamiltonians described in Section III B, we reproduce the tight-binding model of Ref. 4, given by with the matrix elements where ξ = k x /2 and γ = √ 3k y /2 and the lattice constant is set to one.

E. Lattice Hamiltonian of monolayer WTe2
Monolayer WT2 2 was recently discovered to be a twodimensional quantum spin Hall insulator [52][53][54][55] in accordance with previous numerical prediction. 56,57Transport experiments found quantized edge conductivity persisting up to 100K. 54This suggests a much larger band gap compared to devices based on two-dimensional quantum wells. 58It remains an open question whether a simple non-interacting lattice Hamiltonian can reproduce these findings.
We use the restricted set of orbitals in Ref. 59 to construct the spinless tight-binding Hamiltonian.The unit cell contains four sites (labeled A d , A p , B d , B p ) with one orbital on each, and has a symmetry group generated by time reversal, inversion and glide reflection.We use the permutation of the sites under the symmetries and the onsite unitary action (in this case ±1 factors) as input.The model includes hoppings of type A p -A p , B d -B d in neighboring unit cells in the x direction, and B d -A p , A p -B p and A d -B d within the unit cell.We reproduce the Hamiltonian family with 7 free parameters also found in the reference: where This Hamiltonian is identical to the one found previously up to transformations of the Bloch basis.
Extending this analysis to include spin and possible spin-orbit coupling terms, we find that there are 7 additional terms allowed by symmetry in a tight-binding model with the same bonds.The detailed results will be published elsewhere. 60

VI. SUMMARY
Analysis of condensed matter systems is commonly based on single-particle Hamiltonians, the symmetry properties and classification of which are crucial to understanding the physical properties.We discussed the general symmetry structure of single-particle Hamiltonian families, and presented methods to find the full symmetry group of a Hamiltonian, and to generate all Hamiltonians compatible with a given symmetry group.Our methods extend to all continuous and discrete symmetries of single-particle continuum or lattice Hamiltonians.
Although we focused on fermionic systems, the framework we presented is generally applicable whenever the form of the symmetry action and the Hamiltonians is the same, e.g. in the analysis of unconventional superconducting pairing, or even Josephson junction arrays.Our algorithms provide a powerful tool in the ongoing classification of symmetry protected topological phases in a wide variety of physical settings ranging from classical mechanics to circuit QED.The Hamiltonian generator can be extended to search for nonlinear effective field theories and interacting lattice models respecting given symmetries.The symmetry finder may also be further generalized to facilitate more involved symmetry analysis by decomposing group representations.We leave these open questions to future work.
We implemented the algorithms in the Qsymm Python package, making them easily accessible.We demonstrated the usefulness of our approach by applying it to a number of relevant modern research topics including graphene, transition metal dichalcogenides and topological semimetals, resulting in several new insights.

Appendix A: Simultaneous diagonalization
We present an algorithm to simultaneously diagonalize a set of mutually commuting normal matrices.The key property that follows from the commutation is that the matrices share eigensubspaces.By transforming to the diagonalizing basis of one of the matrices, the rest of the matrices are guaranteed to be block diagonal with blocks corresponding to the degenerate eigensubspaces of the first.Considering this, we apply the following recursive algorithm to find the simultaneous eigenvectors spanning the simultaneous eigensubspaces of commuting matrices H i : • If the matrices are 1 × 1, return the 1 × 1 identity matrix.
• Diagonalize the first matrix H 0 , find the orthonormal sets of eigenvectors spanning each (approximately) degenerate eigensubspace.This results in a set of projectors P j onto the eigensubspaces, each consisting of a set of orthonormal column vectors, the number of columns equal to the degeneracy of the j'th eigensubspace.
• If there are no more matrices, return this basis.
• Project the rest of the matrices into each degenerate eigensubspace j: Hij = P † j H i P j for i > 0. • Perform this algorithm on the projected matrices Hij (i > 0) in each eigensubspace j, this returns a set of projectors Pjk .
• Return the set of projectors P jk = P j Pjk for every j and k.
The output is a set of projectors P i , each consisting of a set of orthonormal column vectors spanning a simultaneous eigensubspace of the H i 's.Horizontally stacking the P i 's gives a unitary matrix U such that U † H i U is diagonal for all i.The algorithm is guaranteed to finish, as at each recursion level both the number and the size of the matrices is decreased.The main source of numerical instability is the decision whether to treat two numerically close eigenvalues as degenerate or not.The algorithm is most stable if matrices that have eigenvalues which are either well separated or degenerate to machine precision are first and those which may have accidental near-degeneracies are last.In the physical problems we consider symmetry operators and projectors are of the first kind, while Hamiltonians are of the second.
To solve this problem, we sparsify the matrix representation of a Hamiltonian family by bringing it to reduced row echelon form.In reduced row echelon form, the first nonzero number from the left in a row is always equal to 1, and is located to the right of the first nonzero entry in the row above.Furthermore, the first nonzero entry per row is always the only nonzero entry in its column, and the number of nonzero entries thus minimal.In addition, bringing a matrix to reduced row echelon form preserves its row span.We obtain the reduced row echelon form by performing elementary row operations on the matrix representation of the family.In floating point precision, this generally leads to numerical instability.However, for the applications we consider, this is not a major obstacle, since the matrices we consider are typically small, and usually only contain nonzero elements that are of the order 1, such that the distinction between zero and nonzero entries is unambiguous.
Conserved quantities, which also form a linear space spanned by a set of matrices, suffer from the same ambiguity.We apply the same algorithm of bringing the matrix whose rows are the flattened generators of conserved quantities to reduced row echelon form in order to bring the generator set to a more human readable form.

6 FIG. 1 .
FIG. 1. Pictorial summary of the methods studied in the paper.The symmetry finder and Hamiltonian generator algorithms form a two-way connection between Hamiltonian families and symmetry groups.

FIG. 2 .
FIG. 2. Schematic representation of the onsite symmetry group of a Hamiltonian family (top).The unitary symmetries form a continuous connected Lie group G0.The discrete symmetries can all be combined with any unitary symmetry, forming disconnected components of the symmetry group.For example, T G0 = {T U : U ∈ G0} contains all antiunitary symmetries which are combinations of the canonical time reversal T and some unitary U .Reducing the Hamiltonian family factors out all the unitary symmetries leaving only the identity element e, resulting in a simpler discrete group structure (bottom).
y, z] form a 2s + 1 dimensional spin representation.The generic Hamiltonian again has the form H = 1 (2s+1)×(2s+1) ⊗ H r , because the J i form an irreducible representation of the rotation group.This Hamiltonian, however, is invariant under any unitary transformation of the form U ⊗ 1 with U ∈ U (2s + 1).

FIG. 4 .
FIG. 4. Lattice structures of Kekule-O (left) and Kekule-Y patterns (right).Weak and strong bonds are marked with single and double lines.The high symmetry unit cell is marked in blue.