Group theory of Wannier functions providing the basis for a deeper understanding of magnetism and superconductivity

The paper presents the group theory of best localized and symmetry-adapted Wannier functions in a crystal of any given space group G or magnetic group M. Provided that the calculated band structure of the considered material is given and that the symmetry of the Bloch functions at all the points of symmetry in the Brillouin zone is known, the paper details whether or not the Bloch functions of particular energy bands can be unitarily transformed into best localized Wannier functions symmetry-adapted to the space group G, to the magnetic group M, or to a subgroup of G or M. In this context, the paper considers usual as well as spin-dependent Wannier functions, the latter representing the most general definition of Wannier functions. The presented group theory is a review of the theory published by one of the authors in several former papers and is independent of any physical model of magnetism or superconductivity. However, it is suggested to interpret the special symmetry of the best localized Wannier functions in the framework of a nonadiabatic extension of the Heisenberg model, the nonadiabatic Heisenberg model. On the basis of the symmetry of the Wannier functions, this model of strongly correlated localized electrons makes clear predictions whether or not the system can possess superconducting or magnetic eigenstates.


I. INTRODUCTION
The picture of strongly correlated localized or nearlylocalized electrons is the base of a successful theoretical description of both high-temperature superconductivity and magnetism (see, e.g., [1][2][3] and citations given there). In almost all cases the appertaining localized electron states are represented by atomic orbitals that define, for instance, partially filled s-, d-, or p-bands.
Another option would be to represent the localized electron states by best localized and symmetry-adapted Wannier functions. In contrast to atomic functions, Wannier functions situated at adjacent atoms are orthogonal and, hence, electrons occupying (temporarily) adjacent localized states represented by Wannier functions comply with the Pauli principle. In addition, Wannier function form a complete set of basis functions within the considered narrow, partially filled band. Consequently, Wannier functions contain all the physical information about this energy band.
Wannier functions tend to be ignored by the theory of superconductivity and magnetism because we need a closed complex of energy bands [Definition 2.1] for the construction of best localized Wannier functions. Such closed complexes, however, do not exist in the band structures of the metals where all the bands are connected to each other by band degeneracies.
Fortunately, this problem can be solved in a natural way by constructing Wannier functions with the reduced symmetry of a magnetic group or by constructing spindependent Wannier functions as shall be details in the present paper. In both cases, interfering band degeneracies are sometimes removed in the band structure with the reduced symmetry.
Against the background of the described characteris-tics of the Wannier functions, our following two observations should not be too surprising: (i) Materials possessing a magnetic structure with the magnetic group M also possess a closed, narrow and roughly half-filled complex of energy bands in their band structure whose Bloch functions can be unitarily transformed into best localized Wannier functions that are symmetry-adapted to the magnetic group M . These energy bands form a "magnetic band", see Definition 6.2.
(ii) Both normal and high-temperature superconductors (and only superconductors) possess a closed, narrow and roughly half-filled complex of energy bands in their band structure whose Bloch spinors can be unitarily transformed into best localized spin-dependent Wannier functions that are symmetry-adapted to the (full) space group G of the material. These energy bands form a "superconducting band", see Definition 7.6.
Though these two observations are clear, their theoretical interpretation is initially difficult. This is primarily due to the fact that the models of localized electrons developed so far, as, e.g., the familiar Hubbard model [14], are tailored to atomic orbitals that represent the localized states during an electronic hopping motion. Within modern theoretical concepts, the Wannier functions often are nothing but a complete basis in the space spanned by the Bloch functions. Thus, their symmetry is often believed to do not tell anything about the physics of strongly correlated electrons.
In the light of this background, we suggest to interpret the special symmetries of best localized Wannier functions within the nonadiabatic Heisenberg model [11,15]. This model of strongly correlated localized electrons starts in a consistent way from symmetry-adapted and best localized Wannier functions that represent the localized electron states related to the hopping motion and defines the Hamiltonian H n of the related nonadiabatic system. On the basis of the symmetry of the Wannier functions, the nonadiabatic model makes clear predictions whether or not H n can possess superconducting or magnetic eigenstates [4,9,10,16]. In this context, the nonadiabatic Heisenberg model no longer uses terms like s-, p-, or d-bands, but only speaks of superconducting or magnetic bands.
In particularly interesting cases, the nonadiabatic Heisenberg model predicts that a small distortion of the lattice or a doping is required for the stability of the superconducting or magnetic state. Thus, in undoped LaFeAsO [8] and in BaFe 2 As 2 [9] the antiferromagnetic state must be stabilized by an experimentally well established distortion [17,18], while in YBa 2 Cu 3 O 6 [7] it is stable in the undistorted crystal. Superconductivity in LaFeAsO [8] requires the experimentally confirmed doping [18][19][20][21]. Also the superconducting state in LiFeAs [22] should be accompanied by a small distortion of the lattice which, to our knowledge, is experimentally not yet confirmed. Superconductivity in YBa 2 Cu 3 O 7 [11], MgB 2 [11] as well as in the transition elements [10] (such as in Nb [16]), on the other hand, does not require any distortion or doping.
In the case of (conventional and high-T c [23]) superconductivity, the nonadiabatic Heisenberg model provides a new mechanism of Cooper pair formation which may be described in terms of constraining forces [16] and springmounted Cooper pairs [24].
Any application of the nonadiabatic Heisenberg model starts with a determination of the symmetry of best localized (spin-dependent) Wannier functions related to the band structure of the material under consideration. In the following we shall summarize and update the group theory of Wannier functions as published in former papers and give a detailed description how to determine the symmetry of best localized Wannier functions if they exist in the given band structure. Though we shall also define the two terms "magnetic" and "superconducting" bands which are related to the nonadiabatic Heisenberg model, the presented group theory is independent of any physical model of magnetism or superconductivity. Consider a closed complex of µ energy bands in the band structure of a metal or a semiconductor. Definition 2.1 (closed). A complex of energy bands is called closed if the bands are not connected by degeneracies to bands not belonging to the complex. Definition 2.2 (closed band). In the following a closed complex of µ energy bands is referred to as a single closed band consisting of µ branches.
The metals do not possess closed bands in their band structures. However, closed bands may arise after the activation of a perturbation reducing the symmetry in such a way that interfering band degeneracies are removed. Such a reduction of the symmetry may be caused by a magnetic structure or by a (slight) distortion of the crystal.
Hence, we assume that the Hamiltonian H of a single electron in the considered material consists of a part H G with the unperturbed space group G and a perturbation H H with the space group H, where H is a subgroup of G, In general, the considered closed energy band of µ branches was not closed before the perturbation H H was activated. Assume the Bloch functions ϕ k,q (r) (labeled by the wave vector k and the branch index q) as the solutions of the Schrödinger equation of H to be completely calculated in the first domain of the Brillouin zone. As in Ref. [25], in the rest of the Brillouin zone the Bloch functions shall be determined by the equation where k lies in the first domain, and in the k space outside the Brillouin zone by the equation K denotes a vector of the reciprocal lattice and H 0 stands for the point group of H. Definition 2.4 (symmetry operators). P (a) denotes the symmetry operator assigned to the space group operation a = {α|t α } consisting of a point group operation α and the associated translation t α , acting on a wave function f (r) according to The Bloch functions ϕ k,q (r) of the closed band under observation can be unitarily transformed into Wannier functions centered at the positions R + ρ i , where the functions are "generalized" Bloch functions [25]. The sum in Eq. (2.6) runs over the N vectors k of the first Brillouin zone (BZ), the sum in Eq. (2.7) over the µ branches of the considered band, R denote the vectors of the Bravais lattice, and the coefficients g iq (k) in Eq. (2.7) are the elements of an unitary matrix g(k), (2.8) in order that the Wannier functions are orthonormal, Definition 2.5 (best localized). The Wannier functions are called best localized if the coefficients g iq (k) may be chosen in such a way that the generalized Bloch functions ϕ k,i (r) move -for fixed r -continuously through to whole k space [25].
As it was already shown in Ref. [26], the Bloch functions ϕ k,q (r) as the eigenfunctions of the Hamiltonian H may be chosen in such a way that they vary continuously as functions of k through the first domain and, in particular, approach continuously the boundaries of the first domain. From Eqs. (2.3) and (2.4), however, we cannot conclude that they also cross continuously the boundaries of the domains within the Brillouin zone or at the surface of the Brillouin zone. Fortunately, this problem is solvable by group-theoretical methods [25,27]. Theorem 4.1 shall define the condition for best localized and symmetry-adapted [Definition (2.7)] Wannier functions.

II.2. Symmetry-adapted Wannier functions
In Ref. [25] we demanded that symmetry-adapted Wannier functions satisfy the equation for the elements α of the point group H 0 of H, where the D ji (α) are the elements of the matrices forming a representation D of H 0 which in most cases is reducible. [It should be noted that the sum in Eq. (2.10) runs over w j (r − R − ρ i ) and not over w j (r − R − ρ j ).] Eq. (2.10) defines the symmetry of Wannier functions in general terms, particularly they may be centered at a variety of positions ρ i being different from the positions of the atoms. However, in the context of superconducting and magnetic bands we may restrict ourselves to Wannier functions centered at the positions of the atoms.
Thus, we assume Under these assumptions [15], the Wannier functions may be labeled by the positions of the atoms,

13)
the matrix representatives D(α) of the representation D in Eq. (2.10) have one non-vanishing element D ij (α) with |D ij (α)| = 1 (2.14) in each row and each column, and -Eq. (2.10) may be written in the considerably simplified form and the subscripts i and j denote the number of the atoms at position T and T ′ , respectively.
where {α|t α } ∈ H and R still denotes a lattice vector.
defined by Eq. (2.15) shall be shortly referred to as "the representation defining the Wannier functions" and its matrix representatives D(α) to as "the matrices defining the Wannier functions". In the following Sec IV we shall give a simple condition [Theorem 4.1] for best localized and symmetryadapted Wannier functions yielding the representations of the Bloch functions at all the points k of symmetry in the Brillouin zone. However, in Theorem 4.1 the representations D defining the Wannier functions must be known. Hence, first of all we have to determine in this section all the possible representations that may define the Wannier functions. In this context we assume first that all the atoms are connected by symmetry. This restricting assumption shall be abandoned not until in Sec. III.4.

1)
where R is a lattice vector.
where R denotes a lattice vector. This equation demonstrates that the matrix D(α) has non-vanishing diagonal elements d i (α) if the space group operation a = {α|t α } leaves invariant the position ρ i of the ith atom. These space group operations form a group, namely the group G ρi of the position ρ i . Hence, the non-vanishing diagonal elements d i (α) of the matrices D(α) form a one-dimensional representation d i of the point group G 0ρi of G ρi . The Wannier functions transform according to P (a)w T (r) = d i (α)w T +R (r) for α ∈ G 0ρi (3.5) [cf. Eq. (2.15)] by application of a space group operator P (a) (where R still denotes a vector of the Bravais lattice). From Eq. (2.10) we may derive the equivalent equation (3.6) or, after shifting the origin of the coordinate system into the center of the function w i (r − R − ρ i ), we receive an equation emphasizing the point-group symmetry of the Wannier function at position R + ρ i . In constructing the representation D defining the Wannier functions we cannot arbitrarily chose the onedimensional representations d i of G 0ρi because they must be chosen in such a way that the matrix representatives D(α) form a representation of the point group H 0 , i.e., they must obey the multiplication rule for all the elements α and β in H 0 .
In what follows we assume that all the groups G ρi are normal subgroups of H. In fact, in all the crystal structures we examined in the past, G ρi was a normal subgroup, be it because it was a subgroup of index 2 or be it because it was the intersection of two subgroups of index 2. Both cases are sufficient for a normal subgroup. We believe that in all physically relevant crystal structures G ρi is a normal subgroup of H. If not, the present formalism must be extended for these structures.
When the groups G ρi are normal subgroups of H, each of the groups G ρi contains only complete classes of H, (3.9) We now show that, as a consequence, all the groups G ρi contain the same space group operations.
is an element of G ρi if a ∈ G ρj . Eq. (3.10) even yields all the elements c of G ρi when a runs throw all the elements of G ρj because we may write Eq. (3.10) in the form showing that we may determine from any element c ∈ G ρi an element a ∈ G ρj . On the other hand, Eq. (3.9) shows that c is an element of G ρj , too. When a runs through all the elements of G ρj , then also c runs through all the elements of G ρj . Consequently, all the groups G ρi as well as all the related point groups G 0ρi contain the same elements.
Thus, we may omit the index i and define and respectively, where G ρi and G 0ρi are given by Definition 3.2.  This theorem is necessary, but not sufficient: even if the matrices D(α) can be completely reduced into the irreducible representations of H 0 then they need not form a representation of the point group H 0 [28]. The complete decomposition of a reducible representation is described, e.g., in Refs. [28] and [29], in particular, see Eq. (1.3.18) of Ref. [29]. Theorem 3.1 leads to three important cases: The representation d may be equal to any onedimensional representation of G 0p subduced from a one-dimensional representation of H 0 .
-Case (ii): If all the representations d i are subduced from two-dimensional representations of H 0 , then one half of the representations d i is equal to d A and the other half is equal to d B , where d A and d B are subduced from the same twodimensional representation of H 0 . In special cases, the two representations d A and d B may be equal, see below.
-Case (iii): "Mixed" representations D consisting of both representations d i subduced from one-and two-dimensional representations of H 0 do not exist.
A further case that the representations d i are subduced from three-dimensional representations of H 0 may occur in crystals of high symmetry but is not considered in this paper.
These results (i) -(iii) follow from the very fact that Eq. (2.15) describes an interchange of the Wannier functions at different positions ρ i . Such an interchange, however, does not alter the symmetry of the Wannier functions.

III.3. Sufficient condition for of the representatives D(α) of D
For α ∈ G 0p the matrices D(α) defining the Wannier functions are diagonal, while the remaining matrices D(α) [for α ∈ H 0 − G 0p ] do not possess any diagonal element. Theorem 3.1 only gives information about the diagonal matrices D(α). Hence, this theorem indeed cannot be sufficient because we do not know whether or not the remaining matrices obey the multiplication rule (3.8).
In this section we assume that the matrices D(α) already satisfy Theorem 3.1 and examine the conditions under which they actually form a (generally reducible) representation of H 0 . In doing so, we consider separately the two cases (i) and (ii) of the preceding Sec. III.2.
No further problems arises when case (i) of Sec. III.2 is realized. In this case, Theorem 3.1 is necessary and sufficient. To justify this assertion, we write down explicitly the non-diagonal elements of the matrices D(α).
Let be δ any one-dimensional representation of H 0 subducing the representation d in Eq. (3.13). If we put all the non-vanishing elements of the matrices D(α) equal to the elements δ(α) of δ, then we receive matrices D(α) evidently multiplying as the elements of the representation δ and, consequently, obeying the multiplication rule in Eq. (3.8).

III.3.2. Case (ii) of Sec. III.2
The situation is a little more complicated when case (ii) of Sec. III.2 is realized. Now, the representations d A and d B in Eq. (3.14) may be distributed across the positions ρ i in such a way that the matrices D(α) form a representation of H 0 or do not. Though we always find a special distribution of the d A and d B yielding matrices D(α) actually forming a representation of H 0 , we have to rule out those distributions not leading to a representation of H 0 , because in the following [in Eqs. (4.1), (6.10), and (7.42)] we need the matrices D(α) explicitly.
Let be ∆ [with the matrix representatives ∆(α)] a two-dimensional representation of H 0 subducing the two representations d A and d B of G 0p . The matrix representatives ∆(α) may be determined, e.g., from Table 5.1 of Ref. [29].
As a first step, ∆ must be unitarily transformed (by a matrix Q) in such a way that the matrices ∆(α) are diagonal for α ∈ G 0p , Now consider a certain distribution of the representations d A and d B across the positions ρ i . Then we may determine the elements of the matrices D(α), if they ex-ist, be means of the formula where the ∆ ij (α) denote the elements of ∆(α).
It turns out that in each case the matrices determined by Eq. However, when in Sec. VI or in Sec. VII.3 we will consider magnetic groups, we need all the representatives D(α) of the representation D explicitly. Fortunately, also when the representations d A and d B are equal, Eq. (3.17) is applicable: in this case their exists at least one diagonal matrix representative ∆(γ) of ∆ with vanishing trace and γ / ∈ G 0p . We may define pairs If we find four pairs of positions, we may look for a second matrix representative ∆(γ ′ ) in ∆ with vanishing trace and γ ′ / ∈ G 0p . Then we may repeat the above procedure and receive again four pairs of positions. Now we associate the two representations d A and d B to the positions ρ i under the provision that always positions of the same pair are associated with the same representation d A or d B .
Finally, it should be mentioned that the elements of the non-diagonal matrices D(α) are not fully fixed (as already remarked in Ref. [27]): In Eq. (3.15) we may use the elements δ(α) of any one-dimensional representation δ subducing the representation d. We receive in each case the same diagonal, but different non-diagonal matrices nevertheless satisfying the multiplication rule (3.8). Analogously, in Eq. (3.17) we may determine the matrices D(α) by means of any two-dimensional representation ∆ subducing d A and d B .
In the following Theorem 3.2 we summarize our results in the present Sec. III.3. If not all the atoms at the positions ρ i are connected by symmetry [Definition 3.1], the representation D defining the Wannier functions consists of representatives D(α) which may be written in block-diagonal form where each block comprises the matrix elements D ij (α) belonging to positions connected by symmetry. Otherwise, when the matrices D(α) would not possess blockdiagonal form, Eq. (2.10) would falsely connect atomic positions that are not at all connected by symmetry. As a consequence of the block-diagonal form, the representation D is the direct sum over representations D q related to the individual blocks, We may summarize as follows. The groups of position G p belonging to different blocks may (but need not) be different. However, we assume that the sum in Eq. (3.20) consists only of blocks with coinciding groups of position. If this is not true in special cases, the number µ of the atoms in Eq. (2.7) must be reduced until the groups of position coincide in the sum in Eq. (3.20). Briefly speaking, in such a (probably rare) case atoms of the same sort must be treated like different atoms.

IV. CONDITION FOR BEST LOCALIZED SYMMETRY-ADAPTED WANNIER FUNCTIONS
Remember that we consider a closed energy band of µ branches and let be given a representation D defining the Wannier functions which was determined according to Theorems 3.2 and 3.3. Then we may give a simple condition for best localized symmetry-adapted Wannier functions based on the theory of Wannier functions published in Refs. [25] and [27]. That means, H k is the FINITE group denoted in Ref. [29] by H G k (and listed for all the space groups in Table 5.7 ibidem). Furthermore, let be D k the µ-dimensional representation of H k whose basis functions are the µ Bloch functions ϕ k,q (r) with wave vector k, and χ k (a) (with a ∈ H k ) the character of D k . D k either is irreducible or the direct sum over small irreducible representations of H k .
We may choose the coefficients g iq (k) in Eq. (2.7) in such a way that the Wannier functions are best localized [Definition 2.1] and symmetry-adapted to H [Definition 2.7] if the character χ k (a) of D k satisfies at each point k of symmetry in the first domain of the Brillouin zone the equation The complex numbers d i (α) stand for the elements of the one-dimensional representations d i of G 0p fixing the given µ-dimensional representation D defining the Wanner functions. We add a few comments on Theorem 4.1.
-In Eq. (4.2) we write n i (a) rather than n i (α) because the group G 0p depends on a = {α|t α }.
-The representation D defining the Wannier functions is equivalent to the representation D 0 , i.e., to the representation D k for k = 0, see Eq. (5.10). These two conditions were satisfied in our former papers.
-The irreducible representations of the Bloch functions of the considered band at the points k of symmetry may be determined from the representations D k as follows: , The existence of best localized symmetry-adapted Wannier functions is defined in Satz 4 of Ref. [25]: such Wannier functions exist in a given closed energy band of µ branches if Eqs. (4.28) and (4.17) of Ref. [25] are satisfied. We show in this section that the fundamental Theorem 4.1 complies with these two equations if the Wannier functions meet the assumptions (i) -(iii) in Sec. II.2.
V.1. Equation (4.28) of Ref. [25] As a first step consider Eq. (4.28) of Ref. [25] stating that best localized and symmetry-adapted Wannier functions may exist only if two representations D k ′ ΣR and D k ′ ΣR are equivalent, where we have abbreviated k ′ ΣR by k denoting a point of symmetry lying in the first domain of the Brillouin zone. Consequently, our first task will be to determine the character of D k as well as of D k , The representation D k as defined in Theorem 4.1 is the direct sum of the representations of the Bloch functions of the considered band at point k. The character χ k (a) of the representation D k is simply given by where the matrices D k (a) are the matrix representatives of D k .
The matrix representatives D k (a) of D k are defined in Eq. (4.26) of Ref. [25], is a vector of the reciprocal lattice. Again we have abbreviated k ′ SR by k denoting a point of symmetry. The matrices S(K) as defined in Eq. (4.13) of Ref. [25] are responsible for a continuous transition of the generalized Bloch functions between neighboring Brillouin zones. The matrices D 0 (α) are the matrix representatives of the representation D k for k = 0 as defined in where a = {α|t α } still denotes an element of the space group H. By definition, the matrix M diagonalizes the matrices S(K) which is possible since all the S(K) commute. Thus, the first factor M * S * (K α )M * −1 in Eq. (5.5) is the diagonal matrix where, according to Eq. (2.7) of Ref. [27], also T is a diagonal matrix with Hence, the first factor in Eq. (5.5) may be written as  Strictly speaking, in Ref. [25] we have proven that the condition (5.1) must be satisfied for the points of symmetry lying in the first domain on the surface of the Brillouin zone. Eq. (4.1) demands that in addition the representation D 0 is equivalent to the representation D which is evidently true, see Eq. (5.10).

V.2. Equation (4.17) of Ref. [25]
As a second step we show that Eq. (4.17) of Ref. [25] does not reduce the validity of Theorem 4.1 but this equation is satisfied whenever the assumptions (i) -(iii) in Sec. II.2 are valid. Taking the complex conjugate of Eq. (4.17) of Ref. [25] and transforming this equation with the matrix M * already used in Eq. (5.5), we receive the equation where R j is a lattice vector which may be different in each equation. In fact, these last µ equations (5.16) are satisfied, see Eq. (2.17).

VI. MAGNETIC GROUPS
Assume a magnetic structure to be given in the considered material and let be We demand that the equation is satisfied in addition to Eq.  Again (cf. Sec. II.2), Eq. (6.6) defines the nonvanishing elements of the Matrix N. Hence, also N has one non-vanishing element in each row and each column, As already expressed by Eq. (6.4), we only consider bands of µ branches which are not connected to other bands also after the introduction of the new anti-unitary operation K{γ|τ }. That means that the considered band consists of µ branches as well after as before the introduction of K{γ|τ }. Hence, the matrix N must satisfy the equations In most cases, we may put the non-vanishing elements of N equal to 1. Within the nonadiabatic Heisenberg model, the existence of a narrow, roughly half-filled magnetic band in the band structure of a material is a precondition for the stability of a magnetic structure with the magnetic group M in this material. However, the magnetic group M must be "allowed" in order that the time-inversion symmetry does not interfere with the stability of the magnetic state [9].

VII.1. Definition
Assume the Hamiltonian H of a single electron in the considered material to be given and assume H to consist of a spin-independent part H i and a spin-dependent perturbation H s , Further assume the Bloch spinors ψ k,q,s (r, t) as the exact solutions of the Schrödinger equation Hψ k,q,s (r, t) = E k,q,s ψ k,q,s (r, t) (7.2) to be completely determined in the first domain of the Brillouin zone. Just as the Bloch functions, they are labeled by the wave vector k and the branch index q. In addition, they depend on the spin coordinate t = ± 1 2 and are labeled by the spin quantum number s = ± 1 2 . Consider again a closed energy band of µ branches which, in general, was not closed before the perturbation H s was activated. Now each branch is doubled, that means that it consists of two bands related to the two different spin directions. Just as in Sec. II.1 we assume that the Bloch spinors ψ k,q,s (r, t) are chosen in such a way that they vary continuously through the first domain and approach continuously the boundaries of the first domain. In the rest of the Brillouin zone and in the remaining k space they shall be given again by Eqs. (2.3) and (2.4) [30] where, however, P (a) acts now on both r and t, see Eq. (7.16).
We define "spin-dependent Wannier functions" by replacing the Bloch functions ϕ k,q (r) in Eq. (2.7) by linear combinations of the given Bloch spinors. Thus, Eq. (2.7) becomes and, finally, the spin-dependent Wannier functions my be written as (7.5) Also the spin-dependent Wannier functions depend on t and are labeled by a new quantum number m = ± 1 2 which, in the framework of the nonadiabatic Heisenberg model, is interpreted as the quantum number of the "crystal spin" [15,31,32]. The sum in Eq. (7.4) runs over the µ branches of the given closed energy band, where µ still is equal to the number of the considered atoms in the unit cell.
The matrices which is unitary, (7.8) in order that the spin-dependent Wannier functions are orthonormal, Within the nonadiabatic Heisenberg model we strictly consider the limiting case of vanishing spin-orbit coupling, H s → 0, (7.10) by approximating the Bloch spinors ψ k,q,s (r, t) by means of the spin-independent Bloch functions ϕ k,q (r). In this context, we should distinguish between two kinds of Bloch states ϕ k,q (r) in the considered closed band: where the functions u s (t) denote Pauli's spin functions u s (t) = δ st , (7.12) with the spin quantum number s = ± 1 2 and the spin coordinate t = ± 1 2 . Eq. (7.11) applies to the vast majority of points k in the Brillouin zone.
(ii) If at a special point k the Bloch function ϕ k,q (r) was basis function for a degenerate single-valued representation before the perturbation H s was activated and if this degeneracy is removed by H s , then Eq. (7.11) is unusable for the sole reason that we do not know which of the basis functions of the degenerate representation we should avail in this equation. In fact, in this case the Bloch spinors ψ k,q,s (r, t) are well defined linear combinations of the functions u s (t)ϕ k,q (r) comprising all the basis functions ϕ k,q (r) of the degenerate singlevalued representation (as given, e.g., in Table 6.12 of Ref. [29]). These specific linear combinations are not considered because, at this stage, they are of no importance within the nonadiabatic Heisenberg model.
In the framework of the approximation defined by Eq. (7.11) the two functions ϕ k,q,m (r, t) in Eq. (7.3) (with m = ± 1 2 ) are usual Bloch functions with the spins lying in ±z direction if f ms (q, k) = δ ms . (7.13) Otherwise, if the coefficients f ms (q, k) cannot be chosen independent of k, the spin-dependent Wannier functions cannot be written as a product of a local function with the spin function u s (t) even if the approximation defined by Eq. (7.11) is valid. Consequently, even in the limit of vanishing spin-orbit coupling, the spin-dependent Wannier functions are neither orthonormal in the local space L nor in the spin space S, but in L×S only, see Eq. (7.9). Thus, also in the case H s → 0 spin-dependent Wannier functions clearly differ from usual Wannier functions characterized by The ansatz (7.5) presents the most general definition of Wannier functions. While their localization can be understood only in terms of the exact solutions of the Schrödinger equation (7.2), the limiting case of vanishing spin-orbit coupling characterized by Eq. (7.11) yields fundamental properties of these Wannier functions leading finally to an understanding of the material properties of superconductors [16,32,33].

VII.2. Symmetry-adapted spin-dependent Wannier functions
We demand that symmetry-adapted spin-dependent Wannier functions satisfy in analogy to Eq. (2.15) the equation for a ∈ H since still the assumptions (i) -(iii) of Sec. II.2 are valid. Merely the third assumption (iii) is modified: now the two Wannier functions w T ,+ 1 2 (r, t) and w T ,− 1 2 (r, t) are situated at the same atom and, consequently, we now put where m = ± 1 2 . The vectors T and T ′ are still given by Eqs. (2.13) and (2.16), respectively. The matrices D(α) = [D ij (α)] in Eq. (7.14) are again unitary generalized permutation matrices, and the subscripts i and j denote the number of the atoms at position T and T ′ , respectively, see Definition 2.6.
The operators P (a) now act additionally on the spin coordinate t of a function f (r, t), P (a)f (r, t) = f (α −1 r − α −1 t α , α −1 t), (7.16) where the effect of a point group operation on the spin coordinate t of the spin function u s (t) is given by the equation [29] The matrix denotes the representative of α in the two-dimensional double-valued representation d 1/2 of O(3) as listed, e.g., in Table 6.1 of Ref. [29]. We have to take into consideration that the doublevalued representations of a group g are not really representations of g but of the abstract "double group" g d of order 2|g|, while the single valued representations are representations of both g and g d [29].
Definition 7.1 (double-valued). Though we use the familiar expression "double-valued" representation of a group g, we consider the double-valued representations as ordinary single-valued representations of the related abstract double group g d , denoted by a superscript "d".
Since the index m of the spin-dependent Wannier functions is interpreted as spin quantum number, we demand that the term in Eq. (7.14) describes a rotation or reflection of the crystal spin. Thus, we demand that also the matrices [d mm ′ (α)] are the representatives of the two-dimensional double-valued representation d 1/2 , Consequently, symmetry-adapted spin-dependent Wannier functions are basis functions for the doublevalued representation of H d 0 which is the inner Kronecker product of the single-valued representation D defined by Eq. (7.14) and the double-valued representation d 1/2 . Thus, the 2µdimensional matrix representatives D d (α) of D d may be written as Kronecker products, (7.21) Definition 7.3 (representation defining the spin-dependent Wannier functions). The single-valued representation D of H 0 defined by Eq. (7.14) shall be shortly referred to as "the representation defining the spindependent Wannier functions" and its matrix representatives to as "the matrices defining the spin-dependent Wannier functions". While usual (spin-independent) Wannier functions are basis functions for the representation D defining the Wannier functions, spin-dependent Wannier functions are basis functions for the double-valued representation in Eq. (7.20).
Also the representation D defining the spin-dependent Wannier functions has to meet the conditions given in Sec. III as shall be summarized in Theorem 7.1. The two spin-dependent Wannier func- at the position ρ i are basis functions for the two-dimensional representation of the double group G d 0p related to the point group of position G 0p . The one-dimensional representations d i in Eq. (7.22) fix the (generally reducible) representation D of H 0 defining the spin-dependent Wannier functions [Definition 7.3]. The matrix representatives D(α) of D still are unitary generalized permutation matrices which must be chosen in such a way that they form a representation of H 0 . We again distinguish between the two cases (i) and (ii) defined in Theorem 3.2.
In addition, Theorem 3.3 must be noted.
Theorem 4.1 does not distinguish between usual and spin-dependent Wannier functions but uses only the special representations of the Bloch functions or Bloch spinors, respectively, at the points k of symmetry. Thus, Theorem 4.1 applies to both usual and spin-dependent Wannier functions if in the case of spin-dependent Wannier functions we replace the little groups H k by the double groups H d k . Just as the groups H k , the groups H d k are finite groups in Herrings sense as denoted in Ref. [29] by H G †k and, fortunately, are explicitly given in Table 6.13 ibidem.
When we consider single-valued representations, then the sum on the right-hand side of Eq. (4.1) runs over the µ diagonal elements d i (a) of the matrices D k (a) in Eq. (5.11). When we consider double-valued representations, on the other hand, this sum runs over 2µ diagonal elements d d i,m (a) of the corresponding matrices [where S * (K α ) is given in Eq. (5.8)] because also S d * (K α ) is diagonal and now there are two Wannier func- We need not to solve Eq. (4.1) directly but we may determine the representations D d k complying with Eq. (4.1) in a quicker way. Eq. (7.24) shows that we may write the matrices D The affiliated single-valued band is a closed band that, generally, does not exist in the band structure of the considered material. That means that the Bloch functions ϕ k,q (r) of the closed band under consideration band generally do not form a basis for the representations D aff k even if Eq. (7.11) is valid, see, e.g., the single-valued band affiliated to the superconducting band [Definition 7.6] of niobium as given in Eq. (7.81).
We may summarize the result of this section in Theorem 7.2. Remember that we consider a closed energy band of µ branches and let be given a representation D defining the spin-dependent Wannier functions which was determined according to Theorem 7.1. The band may only be closed after the spin-dependent perturbation H s was activated.
Let be k a point of symmetry in the first domain of the Brillouin zone for the considered material and let be H d k the little double group of k in Herrings sense. That means, H d k is the FINITE group denoted in Ref. [29] by H G †k and explicitly given in Table 6.13 ibidem. Furthermore, let be D d k the 2µ-dimensional representation of H d k whose basis functions are the 2µ Bloch spinors ψ k,q,s (r, t) with wave vector k. D d k either is irreducible or the direct sum over double-valued irreducible representations of H d k . The representations D d k follow Eq. (7.27), where the µ-dimensional representations D aff k define the affiliated single-valued band. Thus, also each D aff k is the direct sum over single-valued irreducible representations of H 0 .
We may choose the coefficients g iq (k) and f ms (q, k) in Eqs. Within the nonadiabatic Heisenberg model we are not interested in spin-dependent Wannier functions that are symmetry-adapted to a general magnetic group as given in Eq. (6.1), but we only demand that they are adapted to the "grey" [29] magnetic group 29) or, in brief, we demand that they are adapted to the timeinversion symmetry. K still denotes the operator of time inversion acting on a function of position f (r) according to Eq. (6.3) and on Pauli's spin functions u s (t) according to Ku s (t) = ±u −s (t) (7.30) (see, e.g., Table 7.15 of Ref. [29]), where we may define the plus to belong to s = + 1 2 and the minus to s = − 1 2 .
The index m of the spin-dependent Wannier functions we still interpret as the quantum number of the crystal spin. Consequently, we demand that K acts on m in the same way as it act on s, Kw T ,m (r, t) = ±w T ,−m (r, t) (7.31) where again we define the plus to belong to m = + 1 2 and the minus to m = − 1 2 . Definition 7.5 (symmetry-adapted to a magnetic group). We call the spin-dependent Wannier functions "symmetry-adapted to the magnetic group M d " as given in Eq. (7.29) if they are symmetry-adapted to H d [Definition 7.2], and if, in addition, Eq. (7.31) is satisfied.
In analogy to Eq. (7.14), Eq. (7.31) may be written as Eq. (7.32) shows that the 2µ-dimensional matrix is the matrix representative of the operator K of time inversion in the co-representation of the magnetic point group derived from the representation D d in Eq. (7.20). Thus, the matrix N d has to comply (Sec. VI) with the three equations (6.9), (6.10) and (6.13) which now may be written as and Within the nonadiabatic Heisenberg model, the existence of a narrow, roughly half-filled superconducting band in the band structure of a material is a precondition for the stability of a superconducting state in this material.

VII.3.2. Time-inversion symmetry of the matrices f (q, k)
In this section we derive the time-inversion symmetry of the matrices f (q, k) defined in Eq. (7.3) and shall give the result in Theorem 7.4. Thought evidence for this important theorem was already provided in Ref. [10] and later papers [16,32], we repeat the proof with the notations used in the present paper.
Within the nonadiabatic Heisenberg model, the validity of this condition is the cause of the formation of symmetrized Cooper pairs in superconducting bands [16,32,33].
This Eq. (7.56) may evidently be written in the more compact form where n is given in Eq. (7.34).
VII.4. k-dependence of the matrices f (q, k) Only those bands are of physical relevance in the theory of superconductivity which are closed not before the spin-dependent perturbation H s is activated. In this section we derive the essential property of such bands and shall give the result in Theorem 7.5.
Let be k a point lying on the surface of the first domain in the Brillouin zone for the space group H and let be H k the little group of k. In this section, k need not be a point of symmetry [according to Definition 4.1] but also may lie in a line or a plane of symmetry. However, we only consider wave vectors k at which Eq. (7.11) is valid. Hence, in general, the Bloch functions ϕ k,q (r) are basis functions for a one-dimensional (single-valued) representation of H k . Nevertheless, in very rare cases, the Bloch function ϕ k,q (r) can be a basis function for a degenerate (single-valued) representation. Both cases shall be examined separately.
Just as in Eq. (3.1) of Ref. [25] we arrange the 2µ Bloch spinors u s (t)ϕ k,q (r) in Eq. (7.11) as column vector with increasing energy, Then the analogous column vector Φ k (r, t) consisting of the Bloch spinors ϕ k,i,m (r, t) in Eq. (7.4) may be written as Φ k (r, t) = g d (k) · f d (k) · Φ k (r, t) The matrices g(k) and f (q, k) are defined by Eqs. (7.6) and (7.7) and still follow Eqs. (2.8) and (7.8), respectively, and We now distinguish between two possibilities: -If the considered energy band was already closed before the spin-dependent perturbation H s was activated, then the affiliated single-valued band actually exists as closed band in the band structure of the material under consideration and, thus, the representations d k,q and d aff k,q are equal, d k,q = d aff k,q . with the consequence that the Wannier functions are, in fact, not spin-dependent but are usual Wannier functions as defined in Eq. (2.6).
-If the considered energy band was not closed before the spin-dependent perturbation H s was activated, then not all the representations d k,q are equal to d aff k,q . Evidently, the qth equation is not solved by f (q, k) ≡ 1 when d k,q = d aff k,q and, consequently, the Wannier function actually are spin-dependent.
We summarize this result in Theorem 7.5.
Theorem 7.5. If the considered energy band was not closed before the spin-dependent perturbation H s was activated, the matrices f (q, k) in Eq. (7.3) cannot be chosen independent of k.
In the Sec. VII.5 the matrix f (q, k) shall by determined for some points in the Brillouin zone of niobium. In rare cases, it can happen that at a special point k some of the Bloch states ϕ k,q (r) are basis functions for a degenerate (single-valued) representation and that this degeneracy is not removed by the perturbation H s . For example, each of the two superconducting bands in the space group P 4/nmm = Γ q D 7 4h (129) listed in Table  3 (b) of Ref. [12] consist of two branches degenerate at points M and A. The single-valued Bands 1 and 2 in Table 3 (a) of Ref. [12] are affiliated to the superconducting Band 1 in Table 3 (b) ibidem; Bands 3 and 4 in Table 3 (a) are affiliated to Band 2 in Table 3 It is crucial for the localization of the spin-dependent Wannier functions that also in this case Eq. (7.64) is solvable. We reveal the solubility of this equation on the example of the bands listed in Table 3 of Ref. [12].
At point M in each of these bands, Eq. (7.69) may be written as Though d kM ⊗ d 1/2 and d aff kM ⊗ d 1/2 again are equivalent, it is not immediately evident that Eq. (7.72) is solvable because f d (k M ) is not a general 4 × 4 matrix. However, also the representations d kM ⊗ d 1/2 and d aff kM ⊗ d 1/2 have a very special form since they may be written simply as Kronecker products. Eq. (7.72) indeed is solvable since it expresses the most general unitary transformation between these special representations.
For instance, consider the point M of one of the bands in Table 3 (b) of Ref. [12] and let be d kM = M 3 given by the calculated band structure of the material under consideration. In addition, let us choose Band 1 in Table  3 (a) of Ref. [12] as affiliated single-valued band. Thus, we have d aff kM = M 2 and Eq. (7.72) is solved by as it may be determined by means of the tables given in Ref. [29]. Though both Band 1 and Band 2 in Table 3 (b) of Ref. [12] are mathematically correct superconducting bands, they cannot be occupied in undoped LaFeAsO [12] which, consequently, is not superconducting.

VII.4.3. Additions
In this subsection we show that neither Eq. (7.51) nor Eq. (7.57) is inconsistent with Eq. (7.64). Remember that in this section we only consider points k at which Eq. (7.11) is valid.
First, taking the complex conjugate of Eq. (7.65), we receive with D aff * k = D aff −k and D * k = D −k the condition showing that we may chose g * (−k) = g(k) (7.76) and, hence, Eq. (7.51) is consistent with Eq. (7.75) and, consequently, with Eq.(7.64). Secondly, transforming the complex conjugate of Eq. (7.69) with the matrix n in Eq. (7.34) and using d * k,q = d −k,q , d aff * k,q = d aff −k,q , and d * 1/2 = n −1 d 1/2 n (7.77) In the same way, we find for the points k F on the line F . Eqs. (7.86) and (7.87) demonstrate that f (1, k) cannot be chosen independent of k in the superconducting band of niobium.

VIII. CONCLUSION
In the present paper we gave the group theory of best localized and symmetry-adapted Wannier functions with the expectation that it will be helpful to determine the symmetry of the Wannier functions in the band structure of any given material. The paper is written in such a way that it should be possible to create a computer program automating the determination of Wannier functions.
In this paper we restricted ourselves to Wannier functions that define magnetic or superconducting bands. That means that we only considered Wannier functions centered at the atomic positions. When other physical phenomena shall be explored, as, e.g., the metallic bound, other Wannier functions may be needed which are centered at other positions, e.g., between the atoms. It should be noted that Refs. [25], [27] and [30] define best localized and symmetry-adapted Wannier functions in general terms which may be centered at a variety of positions ρ i being different from the positions of the atoms.