{'fraction_non_alphanumeric': 0.04815858554663372, 'fraction_numerical': 0.009824737492279184, 'mean_word_length': 4.2369352067118164, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Systematic generation of Hamiltonian families with dualities

Dualities are hidden symmetries that map seemingly unrelated physical systems onto each other. The goal of this work is to systematically construct families of Hamiltonians endowed with a given duality and to provide a universal description of Hamiltonians families near self-dual points. We focus on tight-binding models (also known as coupled-mode theories), which provide an effective description of systems composed of coupled harmonic oscillators across physical domains. We start by considering the general case in which group-theoretical arguments suffice to construct families of Hamiltonians with dualities by combining irreducible representations of the duality operation in parameter space and in operator space. When additional constraints due to system specific features are present, a purely group theoretic approach is no longer sufficient. To overcome this complication, we reformulate the existence of a duality as a minimization problem which is amenable to standard optimization and numerical continuation algorithms. Combined with existing procedures to physically implement coupled-resonator Hamiltonians, our approach enables on demand design of photonic, mechanical, thermal, or electronic metamaterials with dualities.


I. INTRODUCTION
Symmetry groups and their representations are the ingredients that scientists and engineers pour in their cauldron [1][2][3]. Where alchemists sought to make only gold, present day researchers were more ambitious and built metamaterials: synthetic materials with properties nowhere found in nature, such as negative indices of refractions or bulk moduli [4][5][6].
But what is really a symmetry? If we approximate a material as a set of regularly arranged points in space (a crystal), there is an easy answer: symmetries are all transformations of space such as translations, rotations, reflections (in general, isometries) that leave this set of points invariant. These are spatial symmetries (or more generally, space-time symmetries). To describe a material, we also need to attach to each point some internal degrees of freedom. These could represent, for instance, the amplitude of the electric field, the temperature, the displacement of the particles, or the state of a spin at this point. Spatial symmetries transform these internal degrees of freedom in a definite way, depending on whether they are scalars, vectors, pseudo-vectors, etc. In addition, the system might possess internal symmetries that transform the internal degrees of freedom at each point independently (but leave the system globally invariant). As an example, the Ising model is invariant under the inversion of all spins.
In general, symmetries are transformations that leave the system invariant. This leaves the door open for symmetries that are neither spatial symmetries, nor internal symmetries. Can such hidden symmetries exist? If so, are they just coincidences? To address these questions, it will be helpful to consider families of systems depending on continuous parameters, instead of one system at a time. Accordingly, we will consider families of linear operators H(p) smoothly depending on a parameter p. For convenience we will refer to H(p) as a Hamiltonian. Let us however emphasize that (i) this operator might describe systems such as mechanical networks, photonic crystals, diffusive systems, Markov chains, etc. as well as mean-field quantum Hamiltonians, (ii) we will not necessarily assume that H(p) is Hermitian. The Hamiltonian generates the dynamics of the system under consideration, through an equation of motion such as i∂ t ψ = Hψ (or a variant, like ∂ 2 t ψ = Hψ). Here, the vector ψ describes the state of the system, such as the electromagnetic field in a photonic crystal, the displacements of masses in an elastic network, or the wave function in a quantum system. In this context, we focus on transformations acting linearly on this state, namely invertible (usually unitary) matrices U mapping ψ to U ψ. Symmetries are transformations that commute with the Hamiltonian, i.e. such that U HU −1 = H. Symmetries and their breaking play a crucial part in determining the properties of a material, including phase transitions [7], band structures, topological phases [8], and physical responses [2,3]. For instance, piezoelectricity or optical activity can only exist in crystals that do not possess a center of inversion [2,3].
A duality acts in the same way as a symmetry, but it also changes the parameters in a certain way: it maps a system to another system. Formally, a duality is defined by the combination of a transformation U (a matrix acting on the ψ) and a smooth function f defined on the space of parameters p, such as U H(p)U −1 = H(f (p)). For instance, the properties of a system under a magnetic field B (playing the role of the parameter p) are usually different when the magnetic field is inverted to −B. Yet, the behaviors of both systems can often be deduced from one another, as exemplified by the Onsager-Casimir relations in non-equilibrium thermodynamics [9]. In this example, Similar notions of dualities were suggested in systems ranging from self-assembled systems [10] to quantum rotors [11]. Rich consequences arise when a duality is combined with other constraints such as symmetries or conservation laws, ranging from degeneracies in the band structure [12][13][14] and topological states [15,16] to symmetries in the phonon response [17] and constraints on the stiffness tensor of an elastic medium [18].
The goal of this paper is to (i) systematically construct families of Hamiltonians endowed with a given duality and (ii) provide a universal description of Hamiltonians families near self-dual points. We are particularly interested in the consequences of dualities on the properties of spatially extended materials (and metamaterials). To handle such systems, we focus on tight-binding models (also known as coupled-mode theories), which describe a collection of coupled harmonic oscillators ranging from atomic orbitals, mechanical resonators to optical cavities, see Fig. 1. These models provide a somewhat universal effective description of linear systems, provided that one doesn't insist that each individual oscillator should necessarily correspond to an identifiable physical entity. Most crucially from an engineering perspective, systematic methods have been developed to solve the inverse problem of designing realizable systems that correspond to a given tight-binding Hamiltonian [19][20][21][22][23].

II. DUALITIES IN CONTINUOUS FAMILIES OF HAMILTONIANS
Our topic requires some conscious attention to the parameters p, and we denote by P the parameter space in which they live. We call H(p) the Hamiltonian of the system with the parameters p, and ask that the function p → H(p) is smooth.
We say that the family of Hamiltonians H(p) is endowed with a duality when [24] whereÛ is a unitary operator and f is a map from the parameter space P to itself [25]. We emphasize that the duality is determined not only byÛ but also by f . The fixed points of the map f are called self-dual points. At a self-dual point,Ĥ(p) ≡Ĥ(f (p)) and the duality reduces to a symmetry. As a whole, the duality acts on the composite space S = End(H) × P , where H is the Hilbert space of physical states and End(H) is the space of linear maps on H that contains all possible Hamiltonians [26]. The action on this composite space is generated by U = (f, ad U ) (the function ad U is defined by ad U (H) = U HU −1 ). The order of an operation O is the smallest integer n such that O n = Id is the identity (the order is set to ∞ when no such integer exists). As a duality U is composed of two pieces f and U , we can define two separate orders. First, let m be the order of the duality map, i.e. the smallest integer such that f •m (p) = p for all p (here, . Second, let n be the order of the duality operator, such that ad •n U = id [27]. In the following, we assume that m = n, and refer to this integer as the order of the duality. When m = n, the Hamiltonian is either endowed with additional symmetries, or there are redundancies in the parameters (i.e., there are equal Hamiltonians for different values of parameters). For simplicity, we do not consider these situations in this work. Besides, we assume that the order n = m of the duality is finite [28].
In principle, multiple duality operations can be present. Consider two pairs U 1 = (f 1 , ad U1 ) and U 2 = (f 2 , ad U2 ) satisfying Eq. (1). Then, U 1 U 2 = (f 1 • f 2 , ad U1 • ad U2 ) must also satisfy Eq. (1). Hence, we can see the duality operations U a as the representations of an abstract duality group G acting on S. In this paper, we will focus on the situation in which a single duality operation with finite order n is present; then, G is the cyclic group C n .

III. GENERATING DUAL FAMILIES ON DEMAND
A. General strategy Given a duality operatorÛ (acting on a certain space that is also given), we wish to systematically find the families of Hamiltonians p →Ĥ(p) satisfyinĝ Ideally, we would like an algorithm whose input is composed of a suitable description of the space of states (e.g., a crystal C with its Bravais lattice and unit cell and the internal degrees of freedom attached to each site) and of the duality operatorÛ . Its output should consist in some kind of basis whose elements can be combined to express all the families of Hamiltonians satisfying Eq. (2). Methods and tools performing this task have been developed for spatial symmetries. This is the main idea of the theory of invariants of solid-state physics [29][30][31][32][33], which builds upon the abstract theory of invariants in abstract algebra [34]. This approach originated in the study of semiconductors, but it applies to any linear systems, including photonic and phononic crystals [35,36]. Besides, several works have recently been devoted to expanding and automating this procedure using computer algebra systems [37][38][39]. Here, we wish to extend this method to dualities. We will only consider functions f that can smoothly be deformed into isometries by a change of coordinate [40]. Further, we assume that the duality map f acts as a linear isometry on parameter space [41]. Namely, we assume that f (p) = F p in which F is an orthogonal (or unitary) matrix. This assumption is not as restrictive as it may seem: it turns out that in simple cases, it always holds for dualities of finite order [42]. It will allow us to use grouptheoretical methods to describe the general structure of dualities. These will be practically advantageous to obtain explicit parameterized forms of Hamiltonian families with dualities. The simplest application of this general strategy is presented in Box 1, and a concrete example is presented in Box 2. The next paragraphs describe the general case.

Box 1 -dualities in a 1D parameter space
The simplest situation occurs when there is a single self-dual point. This situation can be reduced to the parameter space P = [−1, 1] in which the self-dual point is at p = 0 (in red), with f (p) = −p. p It turns out that the only other possibility is that the entire interval is self-dual (with f (p) = p). Consider on the one hand arbitrary odd and even functions a ± on parameter space, satisfying and on the other hand, self-dual and anti-self-dual matricesĤ ± satisfyinĝ Then, the combinations a + (p)Ĥ + and a − (p)Ĥ − satisfy the duality relation! We can sum any of such combinations, and the most general family of Hamiltonians with duality then readŝ in whichĤ i ± are linearly independent solutions of Eq. (4), and a i ± (p) are arbitrary functions satisfying Eq. (3).
Recall from Sec. II that we have introduced an abstract duality group G. Consider a basis of operatorsĤ µ γ that transform under an irreducible representation (irrep) Γ γ of G (here, µ labels the different basis matrices of the fixed irrep Γ γ ). Hence, in which the matrix ρ(γ) depends only on the irreducible representation Γ γ , and there is an implicit sum over the repeated indices µ . Similarly, consider a basis of functions on parameter space a µ γ : P → P that transform under the conjugate irrep Γ * γ . This means that a µ γ (f (p)) = ρ(γ) µµ a µ γ (p) (8) in which the overline corresponds to complex conjugation. (Note that under our assumptions, the duality map f is a linear map f (p) = F p in which F is a matrix acting on parameters.) We can now combine both ingredients intô We then findÛĤ using that ρ(γ) is unitary [43]. Hence, we have constructed a family of functionsĤ(p) that satisfies the duality equation (1) with f (p) = F p. Of course, one could have multiple sets of basis functions (say, polynomials of different orders), and of basis operators. We have to sum all the combinations with arbitrary coefficients (but the coefficients have, of course, to be the same inside each individual combination).
We can now combine theĤ γ corresponding to all the irreps Γ γ of G intoĤ which is our generic family of Hamiltonians with duality.

B. Cyclic duality group
Up to here, our discussion applies to any finite group G (and with a few modifications, to any compact group).
Action of an order-three duality on parameter space. Here, we consider a duality action f of order three acting on the 2D plane. We have f (2) = id and f (3) = id. In this situation, there is a single self-dual point pc (in red). A typical point p (in blue) has an orbit consisting of three points p, f (p), and f (2) (p). Applying the duality action a third time maps back to the original point f (3) (p) = p.
We now focus on the situation of a single duality operator with finite order n ∈ N. In this case, G is the cyclic group C n ( Fig. 2 shows an example of order n = 3 duality).
The irreps are one-dimensional and the ρ(γ) are complex numbers of modulus one. The irreducible representations of C n are one-dimensional, and have characters Z m n ≡ e im2π/n . Accordingly, the basis matricesĤ m for each irreducible representation are characterized bŷ (Note that this equation could be rewritten without complex numbers, by pairing elements with conjugate characters.) Following the strategy delinated above, we wish to cancel the phase factor e im2π/n by changing the value of the parameters. To do so, consider the two-component parameter p = (p x , p y ) that we represent as the complex number z = p x + ip y . Let us find functions z → a m (z) satisfying a m (e −im2π/n z) = e −im2π/n a m (z).
Then, the combinationŝ satisfy the duality relation Box 2 -self-dual and anti-self-dual Hamiltonians Consider the unitary matrix We wish to obtain the families of Hamiltonians H(p) satisfying U H(−p)U −1 = H(p). To do so, we first seek a basis for the linear spaces of solutions of the (anti-)self-duality equations U HU −1 = ±H. We find two self-dual basis matrices and two anti-self-dual ones, Then, we look for functions satisfying Eq. (3). Let us use polynomials: a ± are respectively polynomials with only even or odd powers of p.
Finally, we combine everything together to write in which all the [a a ± ] i are independent (arbitrary) complex numbers.

C. Spatially periodic metamaterials
We now focus on spatially periodic systems (crystals) which are obtained by repeating a unit cell over a Bravais lattice Γ (see Appendix A for details). To do so, we consider Hamiltonians of the form in whichT (γ) is translation operator by a Bravais lattice vector γ, h(γ) is a finite-dimensional matrix acting on the degrees of freedom in the unit cell [44]. All translation invariant Hamiltonians can be written in this form. For convenience, we also assume that the connections between the oscillators have a finite range (this hypothesis is not crucial, but simplifies the manipulations). Correspondingly, the sum is restricted to a finite subset Γ 0 ⊂ Γ of the Bravais lattice Γ, centered at the origin. We also decompose the duality operator aŝ in which u(γ, γ ) are again finite-dimensional matrices, and in which |γ represents the function in 2 (Γ) having value one at γ and zero elsewhere [45].
Equation (2) can be decomposed over the translation operators and transformed to a finite number of linear relations between the matrix elements h ij (γ, p) and the dual matrix elements h ij (γ, f (p)), gathered in a linear equation of the form (22) in which the matrix elements h ij (γ, p) are all gathered in the vector h(p), etc. and in which A[Û ] and B[Û ] are known matrices (that depend on the duality operator) [46]. Equation (22) can be solved (numerically or symbolically, see e.g. Refs. [47,48]), giving h(f (p)) as a function of h(p) (or conversely).

D. Obtaining explicit parameterized families
The general strategy developed up to here doesn't make any assumptions on the duality, but it leaves us with an set of equations relating h(f (p)) and h(p). In practice, it is convenient to obtain explicit parameterized forms of the Hamiltonian families. To do so, we will restrict our attention to dualities in which the duality map f is linear (as well as all dualities that can be obtained from those by reparameterization). This hypothesis allows us to use the structure described in Section. III A to obtain an systematic way of explicitly constructing the families of Hamiltonians with a duality [49].
To do so, we will first make change of variables to make parameter space as symmetric as possible (we assume that the duality map f becomes an isometry); second, we need to use a basis of functions adapted to the symmetry. We then specialize the strategy of Section. III A to Eq. (22). To get an explicit form, one might further parameterize these functions, e.g. through a series expansion, a Fourier or Chebyshev decomposition, etc. (see an example with polynomials in Boxes 2 and 1).

Box 3 -explicit parameterization for dualities on the line.
We take as a parameter space the real line P = R, and f (p) = −p. There is a single self-dual point at p = 0. Consider on the one hand arbitrary odd and even functions a ± on parameter space, satisfying and on the other hand, self-dual and anti-self-dual vectors h ± , satisfying (Note that h ± do not depend on the parameters.) Then, the combinations satisfy the duality relation. The most general family of Hamiltonians with duality then reads In this equation, h i ± form a basis of independent solutions of the linear Eq. (24), and a i ± (p) are arbitrary functions satisfying Eq. (23) (i.e., odd and even functions). We can then parameterize the functions of a i ± (p), for instance by performing a series expansion near the self-dual point, leading to and (Here, α i n is a coefficient but p 2n is the 2nth power of p.) Let us summarize: we first have to solve Eq. (24) to obtain the basis elements h i ± . Then, we can choose arbitrary coefficients α i n and β i n to get a family of Hamiltonians depending on the control parameter p that satisfies the duality equation. In Sec. III F, we illustrate this strategy on a concrete example of duality.

E. Bloch representation
In spatially periodic systems, it is convenient to use the Bloch (momentum space) representation. We refer to Appendix A for details. In short, each Hamiltonian operatorĤ(p) (for a fixed p) is mapped to a continuous family of matrices k → H(p, k) with the quasi-momentum k in the Brillouin zone (a torus T d ). We emphasize that the quasi-momentum k does not play the same role as the external parameter p: the entire family k → H(k) for k in the Brillouin zone describes a single physical system, while different values of p correspond to different physical systems. In general, the duality changes the quasi-momentum k to Ok, in which O is an orthogonal matrix. We can then write the equivalent of Eq. (1) for the Bloch Hamiltonians [50] F. Example: duality in a 1D crystal Consider first a 1D crystal with two degrees of freedom per unit cell, and the duality defined in Bloch space by f (p) = −p, Ok = −k in Eq. (29), as well as (Here, k ∈ S 1 is the one-dimensional quasi-momentum.) In real space, this corresponds to the operator The second part is due to the action Ok = −k on momentum space. We now consider Bloch Hamiltonians of the form The duality condition Eq. (29) applied to Eq. (32) gives To handle Eq. (33a), let us follow Box 1 and look for a basis H i ± of (anti-)self-dual solutions to U H i ± U −1 = ±H i ± . This has been done in Box 2, and H (0) (p) is given by the RHS of Eq. (19). To satisfy Eq. (33b), we can simply take the matrix elements H Let us be more ambitious and add an additional specification: we now want H(p) to be Hermitian. This requires in addition to Eqs. (33a) and (33b).
Let us first worry about H (0) . We can first impose the constraint Eq. (33a) as explained above. As the generators H a ± in Eqs. (17)(18) are Hermitian, the Hermiticity of H (0) is then ensured by taking the coefficients [a a ± ] i ∈ R in Eq. (19). The second condition requires a bit more work: we are not anymore free to choose an arbitrary function H (1) (p). Instead, we have to find all families of matrices satisfying and then set Proceeding in a similar fashion as in Box 2, we look for a basis B ± = (H j ± ) j of the vector spaces of matrices H ± satisfying (Note that these H j ± are not the same as the ones in Box 2.) where σ j are Pauli matrices. These matrices are then combined with odd/even polynomials (with real coefficients) to construct H (1) (p).
A duality with no degeneracy. We have used a random polynomial of order 2 (random coefficients [a a ± ]i in Eq. (19) and its equivalents). In this example, there ain't no bosonic TRI for to give you no degeneracy. Putting everything together, we obtain Hermitian Hamiltonians Eq. (32) satisfying the duality condition Eq. (29) with the duality operator defined by Eq. (30). In Fig. 3, we show the band structure of one of these Hamiltonians, chosen at random. In this figure, we note that (i) the band structures of dual systems (a) and (c) are obtained by applying the transformation k → −k and (ii) the self-dual band structure is symmetric under k → −k. As we shall see in the next paragraph, richer consequences arise when the duality is combined with another constraints.

G. Other Short Stories
We now show the results of the procedure described in section III A. Before, we must make a few choices: what duality operatorÛ , what parameter space P , and what duality map f shall we impose? We are free to choose [51] any unitaryÛ and any homeomorphism f . The most simple dualities are associated with a one-dimensional The systems (a) and (c) are dual to each other, while the system (b) is self-dual. We consider a 3D cubic lattice with four degrees of freedom per unit cell. Each unit cell is coupled to the six nearest unit cells with an initially arbitrary 4 × 4 hopping matrix. The duality constraint is then applied.
In addition, the coefficients of the Hamiltonian are choosen in order to have a degeneracy at k = 0. The Hamiltonians H(p) in the family are endowed with the duality in Eq. (37) and with bosonic TRI (the physical space Hamiltonians are real-valued). Because of the combination of the self-duality with bosonic TRI, the spectrum (b) is two-fold degenerate.
parameter space P = R and f (p) = −p. We will restrict our attention to this case. We would also like to find duality operatorÛ as simple as possible, but we ask the antiunitary operatorÂ =Û K (obtained by combining the duality operatorÛ with complex conjugation K ) squares toÂ 2 = −1. In this paragraph, we will constrain all tight-binding Hamiltonians to be real-valued (equivalently, when they commute with K ). These requirements guarantee (through Kramers theorem) that the spectra of self-dual Hamiltonians are two-fold degenerate everywhere in the Brillouin zone.
This leads us to consider the duality transformation In momentum space, this corresponds to Eq. (29) with and g(k) = −k.
As we have seen in the previous example (Sec. III F and Fig. 3), it is possible to obtain families of Hamiltonians (with dualities) that have no striking feature: the bands are typically disconnected from each other. When additional symmetries or constraints are imposed (or by accident), it is also possible to have degeneracies (like Dirac or Weyl points) in a self-dual band structure. We illustrate this feature in Fig. 4 for a 1D system, and in Fig. 5 for a 3D system.

IV. DUALITIES IN CONSTRAINED LINEAR SYSTEMS
In the previous sections, we have assumed that all tightbinding Hamiltonians are available to us, perhaps up to some linear constraints such as asking for Hermitian or real-valued Hamiltonians. In particular, we have heavily relied on the hypothesis that linear combinations of physically relevant Hamiltonians are still physically relevant Hamiltonians. Unfortunately, this is not always the case. As an example, consider the dynamical matrixD describing the vibrations of a set of massive particles arranged on a d-dimensional crystal and ruled by Newton equations in which φ are the displacements of the masses with respect to their equilibrium positions (see section IV C below for details). Even though the dynamical matrices D 1 and D 2 might describe two perfectly reasonable networks of particles connected by springs, there is no reason why D 1 + D 2 should necessarily describe yet another networks of spring-connected particles [52]. More generally, one might want to consider systems in which only a certain portion of parameter space is accessible, either for fundamental or for practical reasons. The purpose of this section is to tackle the analysis of dualities in this class of systems. Deprived from the powerful tools of linear algebra and group theory, we are left with no other choice than to resort to numerical analysis. To do so, we start by reformulating the duality condition (1) as an optimization problem. First, we define which is a function defined on the doubled parameter space P × P . Indeed, Eq. (1) is satisfied whenever We define s = (p, p ) and interpret F(s) ≡ F(p, p ) as a vector composed of all its matrix elements [53]. Conversely, any s such that F(s) = 0 is a pair of dual parameters. We can promote this property as a definition, and (re)define self-dual points p as dual pairs of the form s = (p, p). This updated definition purposefully leaves out the function f in Eq. (1). This reformulation affords us two things. First, we can look for pairs of dual parameters using optimization or root-finding algorithms (see e.g. Ref. [48]) by minimizing (or directly finding roots of) the duality potential Second, we can find neighboring pairs of dual parameters from known ones using continuation (path following) algorithms [54,55]. As a particular case, we can also look for self-dual points and continue them by restricting to diagonal parameters s = (p, p).

A. Minimization: finding dual pairs from scratch
We simply look for local minima of the duality potential Importantly, we are only interested in the minima s * that saturate the bound, i.e. those for which F(s * ) = 0. It is possible to use standard minimization algorithms and discard the other minima, or to directly use rootfinding algorithms (see e.g. Ref. [48]). An illustration in a simple mechanical system is presented in Fig. 8 of the next section.

B. Continuation: finding dual pairs from a known one
When a pair of dual points is known, it is possible to use numerical continuation algorithms [54,55] in order to find neighbouring ones. The general strategy is illustrated in Fig. 6. Start from a point s n such that F(s n ) = 0, we perform a series expansion in which δs = s − s n , the index I labels the components of the vector F, the index j labels the components of the doubled parameter s, and we have defined evaluated at s n . As F I (s n ) = 0 by hypothesis, a firstorder prediction of the next point of the curve is s n+1 = s n + δs (47) in which is a small parameter characterizing the size of the steps, and δs is a normalized solution of the linear equations F j I δs j = 0 for all I.
Note that the space of solutions of this equation may have more than one dimension (this occurs if the manifold of solution is not a curve, but e.g. a surface, or when there are branching points). In this case, one has to choose a direction, and different choices will likely lead to different results. The points n+1 (in blue in Fig. 6) is generally not a solution (namely, F(s n+1 ) = 0). To get an actual solution (if it exists), we use a correction step that minimizes F(s) 2 starting from the initial points n+1 (e.g. with Newton iterations). In this way, we obtain a new point on the curve (in red in Fig. 6) This process can then be iterated. Starting from a point s0 in this curve (so F (s0) = 0), we want to find the whole curve. This is the job of numerical continuation (or numerical path following) algorithms. The main idea consists in combining (i) a prediction step (in blue) that starts from a point sn on the curve to suggest a new pointsn+1 that is a bit off the curve, but in the good direction with (ii) a correction step (in red) that projectssn+1 to a new point pn+1 actually on the curve.

General considerations
Consider a set of massive particles arranged on a ddimensional crystal and interacting via short-range pairwise potentials represented by springs. The time evolution of the positions x and momenta π of the particles is described by the (classical) Hamiltonian The Hamiltonian H = T + V is invariant under isometry, but its equilibrium state (x eq i , π eq i = 0) defining the crystal of interest spontaneously breaks these symmetries. (The potential V in H and Eqs. (51-54) is the second term in Eq. (50). It has no relation with the duality potential defined in Eq. (43).) To analyze the linear vibrations of the elastic crystal (phonons), we linearize Hamilton canonical equations (or Newton equations) over the equilibrium state. To do so, let us compute the first derivatives of the potential in which x ij = x i − x j andx = x/ x . These vanish at equilibrium (else there would be a net force on some particle and no equilibrium). The second derivatives are for i = j, and in which r ij = x ij , and for i = j When the springs are at rest in the equilibrium configuration (r eq ij = ij ), the second term in Eq. (52) vanishes and we simply have evaluated at the equilibrium configuration (we omitted the eq labels in the RHS for readability), with the diagonal part is given by Eq. (53), which ensures that the sum rule is satisfied, which is a manifestation of the translation invariance of the original Hamiltonian and ensures the presence of the Nambu-Goldstone modes originating from its spontaneous breaking (acoustic phonons). The displacements u i = x i − x eq i from each particle's equilibrium position then satisfy the equation of motion It is convenient to define φ µ i = √ M i u µ i so that Eq. (56) becomes where we have defined the dynamical matrixD by in which we have assumed that all the masses are positive. Upon choosing a fundamental domain, we can decompose each point x eq i into x eq i = γ + δ n in which γ is a Bravais lattice vector and δ n an element of the fundamental domain. Then in which D(γ, γ ) = D(γ − λ, γ − λ) for all Bravais lattice vectors λ because of the invariance under Bravais lattice translations.
In momentum space (using convention (2) of Appendix A), we get the momentum space dynamical matrix Equations (54) and (55) show that not all matrices can be interpreted as the dynamical matrix of a system of particles connected by springs (even when restricting to symmetric matrices). What's more, the masses and spring constant should be positive in normal circumstances, further restricting the space of allowed matrices.

Example
Consider a one-dimensional mechanical crystal in which the particles are constrained to move on the plane, with two particles per unit cell, as represented in Fig. 7a. Two inequivalent bonds B 1 and B 2 per unit cell connect the nearest neighbors (in purple and green in Fig. 7a). The Bravais lattice is spanned by a single primitive vector a 1 = (1, 0) ∈ R 2 (in red in the figure). To focus on the geometry of the mechanical system, all masses and spring constant are taken to be equal (and set to unity). The mechanical crystal is then fully described by the positions of the two inequivalent particles in the unit cell. We can always set one of the particles at the origin so that its position is (x 1 , y 1 ) = (0, 0). Hence, we are left with the parameters p = (x 2 , y 2 ).
Following the procedure of section IV, we seek self-dual points for the duality operator in whichÎ is defined as in Eq. (38). To do so, we minimize the self-duality potential V = F(p, p) (represented in Fig. 8). We then apply our numerical continuation procedure from one of the self-dual configurations to find a line of dual configurations in the doubled parameter space P × P (i.e., to find pairs of dual points).
The results are presented in Figs. 7 and 8, and can be verified analytically. Expressing (x 2 , y 2 ) as a function of the angle θ between the bonds B 1 and B 2 , we find that the dynamical matrixD(θ) satisfies the duality relation The system is self-dual when φ = 0, i.e. when there is a right angle between the bonds B 1 and B 2 . Note that the duality operator in Eq. (61) is not a spatial symmetry of the system, even at the self-dual point. Typical band structures are plotted in Fig. 7b-d.

V. CONCLUSIONS AND OUTLOOK
In this article, we have developed a theory of dualities in linear systems, with a particular focus on metamaterials. Our results apply to any physical structure described by a linear dynamical system that depends continuously on some parameters.
Using group-theoretical methods, we have shown how to describe generic families of continuous Hamiltonians near self-dual points. These families are generated by generalizations of self-dual Hamiltonians corresponding to the different irreducible representations of an abstract duality group (such as self-dual and anti-self-dual Hamiltonians when the group is Z 2 ). Equipped with this description, we have described a procedure to systematically construct families of Hamiltonians with dualities. Linear constraints on the Hamiltonians such as Hermiticity can be straightforwardly handled in this procedure, which is however limited to the case where non-linear constraint are absent. To handle these, we have reformulated the presence of a duality as a root-finding problem; while this problem is not particularly easy, powerful numerical methods have been developed to tackle it, and can directly be applied.
These procedures apply to any duality, but they don't prescribe what the duality operator should be, nor what should be its action on parameter space. Some practical guidelines to have a two-fold degenerate band structure can be obtained from Wigner's classification of antiunitaries (Appendix C), but this doesn't fully address the issue. This raises the question of whether an enumeration of all possible dualities is possible. In applications to metamaterials, the range of accessible parameters is often limited by practical (not fundamental) constraints. Hence, we have focused on a local approach in parameter space that ignores most of the global structure. In principle, the global topology of parameter space should however affect what dualities are allowed.
While we have focused on the mathematical structures using simple examples, we emphasize that the tightbinding Hamiltonians in our analysis can actually be realized in different kinds of metamaterials. Implementation procedures for tight-binding models recently developed in the context of topological insulators and now routinely used [19][20][21][22][23][56][57][58] will allow to do so, and to explore the consequences of dualities unique to each domain.

Appendices Appendix A: Crystals and Bloch representations
A d-dimensional crystal C is a set of points in Euclidean space that are invariant under a group of translations Γ Z d spanned by the primitive vectors a 1 , . . . , a d of the Bravais lattice. We assume that the states of the system of interest (vibrational states, quantum wave function, etc.) in real space are described by the Hilbert space H of functions ψ : C → V mapping the crystal C to a vector space V that describes internal degrees of freedom at each site. For instance, V could be the vector space R D−1 describing the displacements of a particle in D − 1 dimensions or the vector space C 2 describing the states of a (pseudo)-spin 1/2. To harness the periodicity of the crystal C with respect to the Bravais lattice translations in Γ, we define the quotient C/Γ of equivalence classes of points, and we choose a set of representatives of the equivalence classes, called a fundamental domain or a unit cell F of the crystal (this choice is not unique). The fundamental domain is simply a finite set of points in the crystal that produce the whole crystal (without duplicates) when copied along Bravais lattice translations a i . This allows to uniquely label each point x ∈ C of the crystal as the sum x = γ + ξ of a Bravais lattice vector γ ∈ Γ (itself labeled by an integer in Z d ) and a point in the fundamental domain ξ ∈ F (itself labeled by an integer in {1, . . . , |F|}. We can then represent a physical state ψ in H by all the values ψ(γ, ξ) ∈ V for γ ∈ Γ and ξ ∈ C. Any such state can be decomposed as where ψ y,α are coefficients, δ is the Kronecker symbol, and (v α ) for α = 1, . . . , dim(V) form a basis of V.
In a spatially periodic system, it is convenient to use Bloch (momentum space) representation. To do so, we now identify H with a vector bundle over the Brillouin zone (the family k → H(k) simply corresponds to a section of the endomorphism bundle). There are at least two usual conventions for the Bloch decomposition differing in whether the phase factor attributed to translations is computed from (1) the Bravais lattice translations or (2) the crystal translations, as discussed in Refs. [59][60][61][62][63][64][65][66][67][68]. Informally, they depend on whether one focuses on the quasi-periodic Bloch functions ψ(x) = e ikx u(x) or on the cell-periodic functions u(x). Indeed, the identification of H with a vector bundle over the Brillouin zone requires the choice a basis (more precisely a frame of sections) k → e j,α (k). Here, we consider two possible choices: in convention (1), we use the basis where v α ∈ V is a basis vector for internal degrees of freedom at each site of the crystal, and where j labels the element δ j of a unit cell F. In convention (2), we use where x 0 is an arbitrary origin. In this case, the sum does not depend on the choice of the unit cell, and hence the basis function e (2) j,α doesn't either: j actually labels elements of C/Γ. In the following, we gather the indices j and α in a composite index I = (j, α), and refer to the basis functions as e We are interested in translation invariant operators acting on H, i.e. those that commute with the action of Bravais lattice translations. For concreteness, a Hamil-tonianĤ will be the poster child of such operators in the following, but the discussion of this appendix holds for any linear operator on H, provided that it commutes with translations. It is convenient to use the spatial periodicity to block-diagonalize the Hamiltonian into finitedimensional matrices that can be analyzed more easily. Bloch theorem is the statement that the operatorĤ maps the fiber over k to itself. Hence, we can represent the operatorĤ as a family of matrices k → H(k) (for k in the Brillouin zone). (Without translation invariance, we could still define an object of the form H(k, k ).) Unsurprisingly, this representation depends on the choice of the basis.
We can then define the Bloch Hamiltonian as the family of matrices in which Here, δ j is a point in the fundamental domain F used in the definition of the basis vectors e Temporarily making the dependence on the fundamental domain explicit, we can also write the change of basis matrix between two convention (1) bases with different fundamental domains F and F as Appendix B: Properties of spatial symmetries In this section, we review standard properties of spatial symmetries and of their Bloch representation. We recall that every spatial symmetry (or transformation) can be represented by a unitary matrix that does not depend on momentum, at the possible exception of a phase factor, that disappears in the action by conjugation. (See Refs. [66,70] where it is also proven.)

Spatial symmetries
The space group G of a crystal C is the group of all Euclidean isometries g ∈ Isom(E n ) that satisfy g · C = C and preserve the crystal (we refer to Refs. [71][72][73][74] for more details). Any space group operation g can be decomposed into a rotation R and a translation t as in which we have introduced Seitz notation (R, t) for the space group operation g = (R, t). The inverse of (R, t) is The action (B1) of a spatial transformation g on the Euclidean space effectively defines this operation. We now wish consider the action of g on function defined on the crystal C and taking value in the Hilbert space H of real space crystalline states. This Hilbert space describe the physical degrees of freedom sitting at each point of the crystal: these can be of different nature (scalar, vectors, tensors, spinors, etc.) depending on the physical quantities involved. To take that into account, one has to choose a representation ρ V of the space group G on the internal degrees of freedom V. For instance, a scalar quantity like a temperature will not change at all (ρ V is the identity), while a vector like an elastic displacement should be rotated, etc. A space group operation g ∈ G acts on a function f : C → H as Consider a set of basis functions f α (x) for the Hilbert space H of real space crystalline states. Then where ρ V (g) is a representation of G on V. We will soon drop the · when the context makes it clear.

Bloch representation
A spatial transformation g = (R, t) acts on points in physical space in the way described by Eq. (B1) (this action defines the spatial transformation). Its action on momentum space is defined without ambiguity by requiring the invariance of the scalar product k, r between a momentum k and a position in physical space r. Hence, k → R −1 k when r → Rr (and because R is an orthogonal matrix, we can rewrite this equation using R −1 = R T ). This means that Bloch states with quasimomentum k are mapped to Bloch states with quasimomentum Rk by the symmetry, or in other words that the fiber over k is mapped to the fiber over g −1 k. The action Eq. (B2) on crystalline states can then be represented in Bloch space by a family of matrices of the form in which g· is the action defined by Eq. (B2).
Let us now compute the matrix elements in (B4) from Eqs. (B1-B2) for the space group operation g = (R, t). The results are given by Eqs. (B7) and (B10).
In convention (1), We can uniquely decompose the point gδ j = Rδ j + t in a point δ σg(j) ∈ F plus a lattice vector b j (g) ∈ Γ as This decomposition defines a permutation σ g of the equivalence classes, that does not depend on the unit cell F . However, both δ σg(j) and b j (g) indeed depend on F . We can then write σg(j)α (Rk, x) because the sum over the Bravais lattice can be performed over Rγ or γ + b j (g) instead of γ [75]. Eventually, we find σg(j)α (Rk, x) (B6) and hence g (1) (i,α) (j,β) = e iφ (1) j (g,k) [ρ V (g)] βα δ i σg(j) δ(k , Rk) (B7) in which we have defined the phase factor which depends b j (g) (and hence on j).
In convention (2), we write In the same way as for convention (1), we find σg(j)α (Rk, x) (B9) and hence g (2) (i,α) (j,β) = e iφ (2) This phase factor is independent of α and j, implying that the Bloch representation of g in the second convention only contains an overall k-dependent phase. As expected from Eq. (A5), we find that 3. Which Bloch-space operations cannot be spatial symmetries?
Spatial transformation have the general form Eq. (B3). Conversely, operations that do not have this form cannot be spatial transformation (spatial symmetries).
In momentum space, using the convention (2) for Bloch representations (as defined in Appendix A), a simple necessary condition must be satisfied by all spatial symetries the Bloch representation of any spatial transformation must be of the form where U 0 is a constant matrix (independent of k), as we have shown in Sec. B 2. In other words, U should contain at most an overall momentum-dependent phase e iφ(k) (that cancels out in H → U HU † ). Conversely, matrices that do not satisfy this property cannot represent spatial symmetries. Such a statement cannot directly be made with the convention (1). Let us give an example and a counterexample of spatial symmetry. Glide symmetry in the critical SSH chain. When the two inequivalent atoms A and B are identical (as well as the corresponding bonds), the glide reflection x + 1/2, y is a symmetry. The glide reflection consists in the combination of a mirror reflection σy ≡ x, y and a fractional translation T (a1/2) ≡ x + 1/2, y (in any order).
a. Spatial symmetry. First, consider a onedimensional SSH chain (Fig. 9). It has two inequivalent bonds (i.e., not related by Bravais lattice translations) and atoms. The system has additional symmetries when the inequivalent bonds (atoms) are identical. More precisely, the symmetry increases from p1m1 (less symmetric) to p2mg (more symmetric). (These are frieze groups.) The frieze group p2mg contains a glide reflection x + 1/2, y. The unitary matrix representing this symmetry in convention (1) is in which the size of the unit cell is set to unity (a = 1), and k = k x is the quasi-momentum along the chain. The glide reflection doesn't change k (R = Id with the notations of the previous section). The change of basis matrix defined in Eq. (A6) is in which δ 1 and δ 2 are the positions of the two atoms in the unit cell, and x 0 an arbitrary origin. We can set, for instance, δ 1 = x 0 . Then, the x component of δ 2 − x 0 is the distance between the two inequivalent atoms along the axis, namely half the length of a unit cell. Hence, Then, we find that in convention (2), the glide symmetry is represented by As promised, the non-symmorphic symmetry is represented by a constant matrix multiplied by a momentumdependent phase.
b. Not a spatial symmetry. Let us now seek a counterexample. Let us now write the general form of the 2 × 2 matrix representing a spatial symmetry in convention (1), namely in which we have used that in convention (2), a symmetry is represented by a constant matrix with arbitrary coefficients (unitarity imposes restrictions on the coefficients), up to a momentum-dependent phase φ(k). In 1D, the matrix R in Eq. (B1) can only be R = ±1. The most general W (k) (see Eq. (A6)) is Now, we have (for R = Id) and (for R = −Id) where α = α(k) = k · (δ 1 + δ 2 − 2x 0 ) − φ(k) and β = β(k) = k · (δ 1 − δ 2 ). We now can construct a counterexample: consider It cannot fit in the general form Eq. (B21) because the diagonal elements have different momentum-space dependencies, nor in the general form Eq. (B22) because the off-diagonal elements have different momentum-space dependencies. This counter-example is contrived: it is way easier to determine whether something can be a spatial symmetry from it's representation in convention (2). However, this requires the knowledge of the position of all elements in the crystal, not only of the Bravais lattice.

Appendix C: Kramers degeneracy theorem and antiunitary operators
In this section, we review a theorem of Kramers about the spectrum of linear operators endowed with a particular antiunitary symmetry. We also recall a result of Wigner about the normal form of antiunitary operators, from which all the antiunitary operators to which Kramers theorem applies can be obtained. See Refs. [76][77][78][79] for more details.

Antiunitary operators a. Definition and properties
A map A between two Hilbert spaces is said to be antilinear when A(λφ + µχ) = λA(φ) + µA(χ) for any vectors φ, χ and complex numbers λ, µ (the overbar means complex conjugation). It is customary to omit the parentheses and write Aφ instead of A(φ). An antilinear map A is said to be antiunitary when Aφ, Aχ = φ, χ . Equivalently, Aφ, Aχ = χ, φ (note the reversed order). (Here, ·, · is the Hermitian scalar product that comes with the Hilbert space.) The composition of two antilinear operators is a linear operator, and the composition of two antiunitary operators is a unitary operator.
Every antiunitary operator A can be written as A = U K, where U is a unitary operator and K is complex conjugation in a given basis (it conjugates everything on its right). The inverse of A = U K is then A −1 = U −1 K.
A unitary change of basis V transforms A = U K intõ A =Ũ K whereŨ = V U V −1 . As a consequence, a global phase doesn't affect the square of an antiunitary operator: (e iφ A) 2 = A 2 for φ real.

b. Wigner normal form
It has been shown by Wigner [76,77] that any antiunitary operator (on a finite-dimensional Hilbert space) A can be written in some basis as a block matrix (i.e., decomposed into the direct sum) A = e iφ1 a 1 ⊕ · · · ⊕ e iφn a n = where α = 0 is real [and depends on the block i] (the case α = 0 can be reduced to A 0 ⊕ A 0 ). These blocks square respectively to A 2 1 = 1 and A(α) 2 = diag(e −iα , e iα ). In particular, A(π) 2 = −Id (Id is the identity matrix).
This decomposition shows that an antiunitary operator A with A 2 = −Id (to satisfy the hypotheses of Kramers theorem, see next section) should be composed of blocks A(π) ≡ iσ y (σ y is a Pauli matrix), and in particular should be even-dimensional.

Kramers theorem
Let us consider an antiunitary operator A such that A 2 = −Id. Let us assume that some linear operator H commutes with A, namely AHA −1 = H. Then, for every eigenvector |ψ of H with eigenvalue λ, the vector A |ψ is an eigenvector of H with eigenvalue λ, and it is orthogonal to |ψ . When H is Hermitian, all its eigenvalues are real (so λ = λ). As a consequence, the linearly independent vectors ψ and Aψ are degenerate. This result is known as Kramers theorem [78,79].
The theorem actually contains two sub-statements. First, every vector ψ is orthogonal to Aψ, because ψ, Aψ = A 2 ψ, Aψ = − ψ, Aψ (C3) using A 2 = −1. This property crucially requires A 2 = −1, but H is out of the picture. Second, consider an eigenvector ψ of H with eigenvalue λ. Then, applying A to Hψ = λψ (C4) leads to AHψ = Aλψ (C5) and using that A and H commute, Hence, Aψ is an eigenvector of H with eigenvalue λ. This statement requires [A, H] = 0, but the value of A 2 doesn't matter.
system; we still have to specify how they act. To do so, we consider an action f : G × P → P of G on the parameter space. Each g ∈ G gives rise to a function p → fg(p) on P . This specifies how parameters change under the duality g. Then, we define a linear action U : G × E → E of G on the vector bundle E of physical states. It specifies how states are transformed. We impose that Ug = U (g, ·) is a linear map from the fiber Ep over p to the fiber E fg (p) over fg(p). This is the equivariance condition, which simply means that a state ψ of the system with parameters p is mapped to a state Ugψ of the system with parameters fg(p). We can deduce how operators transform from the way states transform. Accordingly, an operator H acting on the states of the system with parameters p is mapped to UgHU −1 g , which is an operator acting on the states the system with parameters fg(p).
[27] We can define a separate order n such that U n = Id is the identity operator. In general, n = n .