A Guide to the Bott Index and Localizer Index
Terry A. Loring
math-ph cond-mat.mes-hall math.MP math.OA






Department of Mathematics and Statistics, University of New Mexico,
Albuquerque, New Mexico 87131, USA

loring@math.unm.edu



The Bott index is inherently global. The pseudospectal index is inherently
local, and so now the preferred name is the localizer index. We look
at these on a rather standard model for a Chern insulator, with an
emphasis how to program these effectively. We also discuss how to
tune the localizer index so it behaves like a global index.


A Guide to the Bott Index and Localizer Index
    Terry A. Loring
    
=============================================





§ INTRODUCTION


The most elegant real-space topological invariant to describe a Chern
insulator is computed on an infinite model, the index of a Fredholm
operator:

    ind(PX+iY/|X+iY|P+I-P).

Here X and Y represent position operators and P is the Fermi
projection derived via functional calculus from the Hamiltonian H.
This does not require translation invariance, just a gap in the spectrum
of H containing the Fermi level. By the work of Bellissard, van
Elst, and Schulz-Baldes <cit.> we know
this equals the Chern number in the case of translation invariant
systems.

For numerical studies we must work with finite matrices. This can
be done in several ways. Some infinite area formulas can be approximated
by finite area calculations. For example, Prodan explains in <cit.>
how a finite area approximation leads to a number that, while not
an integer, converges rapidly to the noncommutative Chern number.
Kitaev <cit.> defined a generalized notion
of the Chern number with a formula to be computed over an infinite
area. Several authors have noted that in many instances a finite-area
approximation converges well to Kitaev's generalized Chern number
<cit.>. This method has the
advantage of giving a local index. Thus two lattice models can be
melded on a linear interface and a local (non-integral) real-space
marker can be defined that seems to sensibly mark topological and
trivial areas <cit.>.

An alternate approach is to define, via K-theory, an integer associated
to a finite model, and attempt to prove that this integer coincides
with the infinite-area index in (<ref>). The
C^*-algebra in this case is typically very elementary, just the
n-by-n complex matrices M_n(ℂ). In
proving things about such an index, one tends to need complicated
C^*-algebras defined by generators and relations <cit.>
and tools like E-theory <cit.>.
Keeping the definitions simple should, in general, lead to fast algorithms.

The essence of this approach is to look to the old theory of what
K-theory can tell us about almost commuting matrices <cit.>.
For example of what was known by the 1990s, consider two unitary matrices
that almost commute. These are sometimes small perturbations of commuting
unitary matrices, sometimes not. It can be shown <cit.>
that being close or far from commuting unitary matrices can be inferred
from the value of a K-theoretic invariant. There are ways to compute
this invariant that do not really look like K-theory <cit.>,
and other formulas that are straight-forward computations involving
formulas for idempotents. Indeed, there are many, many formulas we
might choose <cit.>. 

Finite models force upon us a choice of boundary conditions. Actually,
first one selects a shape, with square being the de facto choice.
A frequent choice in physics is to use periodic boundary conditions,
which implicitly turns the square into a torus. This is great mathematically,
being an enabler of the Fourier transform. However, this choice eliminates
edge modes, arguably the most attractive feature of a Chern insulator.
The other choice is working with open boundaries, where on a lattice
model the hopping terms in the infinite-area Hamiltonian are dropped
if they involve sites outside the finite patch. Mathematically, this
is just the obvious compression of the infinite-area Hamiltonian.
This leads to edge modes in the finite model, and is more realistic
than a torus model. This is actually an annoyance in some instances,
as the edge modes can make it harder to see bulk phenomenon in a small
system <cit.>.



§ THE CLIFFORD SPECTRUM


The Clifford spectrum gives a nice way to see the difference between
the Bott index and the localizer index. In both cases there is a topological
space associated to a finite system, but the spaces are very different.
All sorts of spaces can emerge as the Clifford spectrum of a finite
collection of Hermitian matrices. This is well established in string
theory, for example in <cit.>.
Kisil <cit.> observed earlier that the joint spectrum
of a few Hermitian matrices could be a surface.

For a periodic system, square with sides of length L, one cannot
easily work with the position operators X and Y as these will
have commutator norm with the Hamiltonian H on the order of L.
Thus one considers “periodic observables”

    U=e^2π i/LX,   V=e^2π i/LY

which really correspond to four observables (Hermtian matrices)

    M_1=cos(2π i/LX),   M_2=sin(2π i/LX),   M_3=cos(2π i/LY),   M_5=sin(2π i/LY).

One has then a fifth matrix 

    M_5=P

where P is the spectral projection for H associated to energy
levels below the Fermi gap.

The Clifford spectrum of Hermitian matrices M_1,….M_d (all
of the same size) is a subset Λ(M_1,…,M_d) of ℝ^d,
defined as follows. One chooses nontrivial Γ_1,…,Γ_d
that form a Clifford representation, here meaning 

    Γ_j^†=Γ_j,  Γ_j^2=I

for all j and 

    Γ_jΓ_k=-Γ_kΓ_j

 whenever j≠ k. Then one forms the localizer 

    L_(M_1,…,M_d)=∑(M_j-λ_j)⊗Γ_j

(called sometimes the localized Dirac operator <cit.>).
Then

    λ∈Λ(M_1,…,M_d) L_(M_1,…,M_d) is singular

 defines the Clifford spectrum of (M_1,….M_d)
, and this is independent of the choice of Γ_j.

For reasonable H (one needs some locality, boundedness and a reasonable
gap) when L is large enough, we have U, V and P that satisfy
the following:




[c]0.7

    U^†U≈ I,   UU^†≈ I,   V^†V≈ I,   VV^†≈ I,


    UV=VU,   UP≈ PU,   VP≈ PV,


    P^2=P,   P^†=P.






All this can be made ever so rigorous, as in <cit.>.
What is critical is that if these relations were exact, they would
become equations in ℝ^5 that define a copy of the space

    𝕋⊔𝕋,

meaning two disjoint copies of the torus. 

What we have then is a fuzzy version of two copies of the torus, and
the Clifford spectrum will be a subset of ℝ^5 that is
close to the union of two tori. It is the topology of 𝕋⊔𝕋
that is used in defining the Bott index. Perhaps, when the Bott index
is nonzero, and in a clean system on a regular lattice, the Clifford
spectrum of M_1,…,M_5 is homeomorphic to 𝕋⊔𝕋.
Proving this is out of reach, for now.

The localizer index is based on different topology, and it starts
with open boundaries. In this case, H will almost commute with
X_0=2LX and Y_0=2LY. Notice we
cannot expect a gap in the spectrum of H. One can spectrally flatten
the infinite-area Hamiltonian and compress that, but the result will
still be ungapped due to edge modes.

Let us assume the system is centered at 0. The Clifford spectrum
of (X_0,Y_0,H) will be some compact subset of ℝ^3.
The key relation here is

    L_0(X_0,Y_0,P)≥ c,

where c>0 is roughly that spectral gap in the infinite area Hamiltonian.
See <cit.> for details on what we know about c.
We also have relations




[c]0.7

    -1≤ X_0≤1,  -1≤ Y_0≤1,


    X_0Y_0=Y_0X_0,   X_0H≈ HX_0,   Y_0H≈ HY_0,


    -C≤ H≤ C






where C is the norm of the infinite-area Hamiltonian. We have here
a fuzzy version of 

    S(C,c)=([-1,1]×[-1,1]×[-C,C])∖ B(c)

where B(c) is the open ball of radius c at the origin. The pseudospectrum
Λ(X_0,Y_0,H) is expected to be some manner of a surface
within S(C,c). While S(C,c) is homotopic to a sphere, the space
Λ(X_0,Y_0,H) will be more complicated, and is very difficult
to compute. Still, it will have a K-theory class coming ultimately
from the K-theory of the sphere which will lead to the localizer
index.

If one uses the spectrally flattened infinite-area Hamiltonian and
compresses that to the finite square system, one has then a Hamiltonian
H_1 that has “spectrum localized near (0,0) equal to a fuzzy
{0,1}” in the sense put forward in <cit.>. A
naive expectation is that the Clifford spectrum of (X_0,Y_0,H_1)
will be roughly two copies of the square. In the case of trivial insulators,
this is what one should to find. It is the edge states, when the system
is topological, that work to glue together the top and bottom square
to create more of a cube. This is all heurisitic, until we find betters
ways to compute the Clifford spectrum. In the case of the Bott index,
the periodic boundary conditions turn the top and bottom square each
into a torus. In the case of the localizer index, when the system
is topologically non-trivial, the edge states add to the Clifford
spectrum, forming walls that join the top and bottom squares into
something like a sphere.



We can, with lots of computer time, compute the Clifford pseudospectrum
of (X_0,Y_0,H). This is essentially a blurred out version of
the Clifford spectrum, and formally is the scalar-valued map on 3-space

    λ↦‖ L_(X_0,Y_0,H)^-1‖ ^-1.

Where this is zero is the Clifford spectrum. We can more easily compute
where this is small, and get an estimate on where lives the Clifford
spectrum. 

Now that we are working with just three observables, we can specify
a choice for the Γ_j, specifically Γ_1=σ_x,
Γ_2=σ_y, Γ_3=σ_z. Then

    L_(M_1,M_2,M_3)=(M_1-λ_1)⊗σ_x+(M_2-λ_2)⊗σ_y+(M_3-λ_3)⊗σ_z

so 

    L_(X_0,Y_0,H)=[[                H-λ_3 (X_0-λ_1)-i(Y_0-λ_2); (X_0-λ_1)+i(Y_0-λ_2)               -H+λ_3 ]].


Ideally, the reader would see this computed for the Hamiltonian considered
in the rest of the paper. As the author has consumed already considerable
resources from the Center for Advanced Research Computing at the university
of New Mexico, it seems prudent to offer the reader a different example,
which was computed earlier. Here the Hilbert space is built from a
square sample of a quasilattice, and the Hamiltonian is a variation
on a “p_x+ip_y” tight binding model, as in <cit.>.
A portion of the pseodospectrum is indicated in Figure <ref>.
The Hamiltonian is the compressed original Hamiltonian, not a compression
of the spectrally flattened Hamiltonian. The pseuodospectrum extends
up to about λ_3=6 and down to about λ_3=-6.

The localizer can have a very large gap even when the Hamiltonian
does not. This happens because, by design, the localizer can focus
on “bulk states” while ignoring edge states. This can happen even
with strong disorder, as is illustrated by the disorder-averaged study
in Section <ref> and by the examination
of the localizer spectrum in <cit.>.



Unlike the Bott index, the localizer index is flexible in what geometry
is used. For example, one can use a square with a disk cut out of
the middle. In this case, there are additional edge modes and the
result is more like a torus than a sphere. See Figure <ref>.



§ FORMULAS FOR THE BOTT INDEX


The first step in finding an invariant for matrices with the relations
as in Equation <ref> is to eliminate the second
torus (the one at high energy). This is the step of considering band-compressed
angular position operators

    PUP:P(ℂ^n)→ P(ℂ^n)

and

    PVP:P(ℂ^n)→ P(ℂ^n).

A more precise notation is .(PUP)|_Pℂ^n
and .(PVP)|_Pℂ^n. These are operators,
not matrices. To obtain matrices we must select a basis for the range
P(ℂ^n) of P. Assuming one obtains the Fermi projector
P by doing a full diagonalization of the Hamiltonian H∈M_n(ℂ)
such a basis is then available at no further cost.

This is a sticking point. Taking a full diagonalization of H forces
us to work with dense matrices. As we are dealing with approximate
relations, we could work with an approximation P_0 to P that
is sparse. This is a standard procedure in quantum chemistry <cit.>.
However, there are other parts to the Bott index formula where a sparse
approach is not at all clear. Also, we need sparse approximate spectral
flattening given only a mobility gap, which will be harder than the
case of a true gap. For now, the Bott index is a dense matrix formula,
which is its major limitation.

Suppose then that you have used an eigensolver on H to find for
ℂ^n a basis q_1,…,q_n
of eigenvalues, with q_1,…,q_m
in the subspace below the Fermi level and q_m+1,…,q_n
above. (In a disordered system, there may be no apparent gap, but
one still is selecting a Fermi level.) Then q_1,…,q_m
is a basis for P(ℂ^n). The matrices, with respect to
this basis, for the operators .(PUP)|_Pℂ^n
and .(PVP)|_Pℂ^n are

    U_1=W^†UW,   V_1=W^†VW,

where

    W=[[            ;            ; q_1   ⋯ q_m;            ;             ]].

Notice W is a partial isometry. 

The smaller matrices U_1 and V_1 satisfy the more familiar
relations




[c]0.7

    U_1^†U_1≈ I,   U_1U_1^†≈ I,   V_1^†V_1≈ I,   V_1V_1^†≈ I,


    U_1V_1≈ V_1U_1






There is a theory on how to create a K-theoretic index for such
approximate relations <cit.>.
Such an index is never unique. All we can say is that for very small
commutators two such indices will agree. We can't usually say how
small the commutator must be for two such invariants to agree. Some
of these choices for invariants will make more sense physically than
others. Some will be faster to compute than others.

The general method here starts with a formula for a projector-valued
function on the torus, and 2-by-2 matrices suffice. For example 

    (w,z)↦[[        f(z) g(z)+ih(z)w; g(z)-ih(z)w      1-f(z) ]]

for |v|=|w|=1. In fact, this need not take values only in projectors:
approximate projectors (p^2≈ p) will be fine. Thus we assume
fg≈0 and f^2+g^2+h^2≈ f. The formulas used
must make sense when applied to noncommuting, nonunitary matrices.
For example, one can take f, g and h to be Laurant polynomials
in z (so trig polynomials in the argument of z). Then the K-theoretic
index is the number of eigenvalues above 12, minus the
number below, of the matrix 

    [[              f(V_1)   g(V_1)+ih(V_1)U_1; g(V_1)-iU_1^†h(V_1)            1-f(V_1) ]].

This is a simplified version of what was called the trig method in
<cit.>.

There are many variation possible on this. Moreover one can work directly
with U, V and a sparse approximation P and use the theory
of <cit.> for the space
𝕋⊔𝕋. This might lead to faster algorithm
than is a available to the Bott index, in theory. 

An issue with all invariants using this method is that they do not
treat all parts of a square lattice equally. This is because the function
in (<ref>) is a continuous function from the
torus to the Bloch sphere. Every such map has at least one point on
the sphere where some nonzero area of the torus is smashed together
sent to a region of zero area on the sphere.

In the special case of the approximate relations in Equation <ref>,
there is an elegant alternative to produce indices proven equivalent
(equal for small commutators) to an index as above. For the case of
almost commuting unitary matrices, this was proven a while back <cit.>.
The formula was modified to work for almost unitary matrices that
almost commuting more recently <cit.>. It also treats
U_1 and V_1 symmetrically, and treats all areas in the lattice
equally (if the lattice is treated as on a torus). This alternate
formula is what physicists now call the Bott index.

The Bott index uses the smaller matrices from Equation <ref>
and is

    Bott(U,V,P)=(1/2π iTr(log(U_1V_1U_1^†V_1^†)))

where  indicates real part. It is guaranteed, in theory, to
makes sense, and be in integer, when (∞,0] is not in the spectrum
of U_1V_1U_1^†V_1^†. In practice eigenvalues
of U_1V_1U_1^†V_1^† very near (∞,0]
also cause trouble. This formula uses the matrix logarithm, and you
don't want to compute a matrix logarithm as that takes a long time.

If X is a matrix with no spectrum lying in the set (-∞,0),
then

    (1/2π iTr(log(X)))=∑(1/2π ilog(λ/|λ|))

where λ_1,…,λ_n are the eigenvalues of X
listed according to multiplicity. In all instances the usual branch
of log is assumed.



This follows easily from the spectral mapping theorem for analytic
functions.

Notice we have no need for any eigenvectors of DCD^†C^†,
just the list of eigenvalues. After multiplying matrices to form the
commutator, you the calculate all the eigenvalues, take the scalar
logarithm of each, and add those. The heart of the algorithm is as
follows, when implemented in Matlab. 











Here  is the partial matrix of eigenvectors and 
and  are the normalize exponentiation position matrices
U and V. Also and  and  are what have been
denoted U_1 and V_1. See the file 
in the supplementary files <cit.>.

It could be that someone has computed the projector P onto the
subspace of states below the gap, but not a orthonormal basis for
that space. This will happen if you use the methods of quantum chemistry
<cit.> to approximately compute the Fermi projector
P. In that case, one can use the almost commuting matrices 

    A=PUP+(I-P) and B=PVP+(I-P)

and it can be shown that

    Bott(U,V,P)=(1/2π iTr(log(BAB^†A^†))).

Since A and B are larger than U_1 and V_1 then this
may end up as the slower method, but perhaps someone knows how to
compute the trace of the log of a sparse matrix quickly. Without a
method to deal with that trace of a matrix log, one is advised to
work only with U_1 and V_1.



§ EXTENDING AN OLD STUDY 




For the purposes of evaluating the algorithm of the Bott index, we
replicate and extend the numerical study in <cit.> (see
also <cit.>). The reason to select a known model
Hamiltonian and disorder-induced transition is to keep the focus here
on the methods for computing the Bott index.

We start with a regular square lattice, with periodic boundary conditions,
with two basis elements at each site. Upon this we consider a real-space
tight binding model, associating a Hamiltonian term to each site of
the tiling and a hopping matrix to each link between closest neighboring
sites. The on-site element corresponding to site j is 

    H_j=(M-4B+W_j)I_2

where M and B are constants, and W_j varies by site, each
W_j independently taken from a uniform distribution in [-4,4].
The full hopping term between sites j and k, for j to the
left of k is

    H_jk=Bσ_z+Aσ_x

for j below k is

    H_jk=Bσ_z+Aσ_y

The constants we will use are

    A=1, B=-1, C=0, M=-2

which is as was done in producing <cit.>. 

In other parts of <cit.> M=-1 was sometimes used.
The meaning here of the constant A is off by 12 how
how it is used in <cit.>. We continue to use the
convention from <cit.> as listed above, and in any case
the reader is encouraged to look at the file in the supplementary files <cit.> to see how the Hamiltonian
and lattice are constructed.



As in <cit.> we vary the Fermi level for many Hamiltonians
with different disorder terms. The index averaged over disorder that
resulted is shown in Figure <ref>. The improvments
in computing hardware over the past twelve years gave a moderate improvement,
as now it is feasable to work up to L=80 instead of only L=60,
as in <cit.>. Also, more samples were computed,
so the resulting plots are clearer.

As we are computing all eigenvalues of some dense matrices, the time
needed to compute the Bott index is expected to grow a the cube of
system size, so 𝒪(L^6). Multiple values of the Bott
index, for multiple Fermi levels, can be computed from the same diagonalization
of the Hamiltonian, so in timing this algorithm it makes sense to
time the Hamiltonian diagonalization separately from the time to compute
a single Bott index given that diagonalization. As one can see from
the left panel of Figure <ref>, both tasks do
take time that grows in what looks like 𝒪(L^6). 



As a comparison, the timing is reported (right panel of Figure <ref>)
as well for the algorithm that uses the larger matrices of Equation <ref>.
There is accumulated numerical error that keeps the result from being
exactly an integer. This is generally ignored; one just rounds to
the nearest integer. In Figure <ref> is shown
the distance to the nearest integer. This seems like it will be negligible
even for systems so large that computing a single Bott index will
be impractical.

It is worth pointing out that we do not, it seems, have a theory that
tells us what the correct answer is here for the disordered averaged
Bott index. However, the Bott index is known to correlate, in the
gapped case, to the Kubo conductance <cit.>.
It has been used and evaluated in many settings, such as <cit.>.
What Figure <ref> is showing is that the computed
Bott index come out very close to an integer, the error coming from
the cumulative errors in calculating eigenvalues numerically.

The plot of the time needed makes it seem that computing a Bott index
for a 100-by-100 system is reasonable. If computing an average over
disorder, and at multiple Fermi levels, this is not the case, as the
total time needed is excessive. On the other hand, since the number
of samples needed decreases with system size, the growth in that case
is more like 𝒪(L^5). That is still a problem, so if
one wants to push this study further, one needs to be doing sparse
matrix computations.



§ THE LOCALIZER INDEX 


The localizer index works only with open boundaries. Unlike the Bott
index, it can easily be modified for different symmetry classes and
physical dimensions <cit.>. Here we focus on
2D system in class AI. Several papers <cit.> referred
to this index as the pseudospectrum index, but that is misleading
as the index is defined where the pseudospectrum isn't.

We have three observables (X,Y,H) with X and Y position and
H the Hamiltonian. We could change units in X and Y and would
still have three observables, now (κ X,κ Y,H) for some
nonzero, finite scaling constant κ. If λ
is not in Λ(κ X,κ Y,H) then the localizer
index at λ for (κ X,κ Y,H)
is

    1/2Sig(L_λ(κ X,κ Y,H)),

where L_λ is the localizer of Equation <ref>
and Sig(M) refers to the number of positive eigenvalues
minus the number of negative eigenvalues of the nonsingular Hermitian
matrix M.

In <cit.> the localizer was referred to as the
Bott element since it does give a representative of the Bott element
when applied to standard position coordinates on the sphere. Now that
the Bott index is rather firmly entrenced in physics, it seems prudent
to call it the localizer. This reflects the fact that one can use
the localizer to differentiate edge from bulk states <cit.>.

It is important to not compute all the positive or all the negative
eigenvalues of the localizer when computing the index. That would
defeat the advantage we gain from the fact that L_λ(κ X,κ Y,H)
should be a sparse matrix. One might need to zero-out small hopping
terms between distant sites to ensure H is sparse. In general κ X
are κ Y are diagonal, and so very sparse.

If M is nonsingular and Hermitian, one computes an LDLT decomposition
of M and uses Sylevester's law of inertial to compute the signature
of M. This means, in practice, on finds a factorization

    M=S^-1PLDL^†P^†S^-1

where S is diagonal with positive values on the diagonal, P
is a permutation matrix, L is lower triangular, and D is block-diagonal,
with both 1-by-1 and 2-by-2 blocks. Since

    M=(S^-1PL)D(S^-1PL)^†

we have 

    Sig(M)=Sig(D)

and the signature of D is the sum of the signatures of the blocks
of D. 



Unfortunately, we require LDLT for a sparse complex Hermitian matrices,
while Matlab <cit.> and MUMPS <cit.> only
allow for sparse real symmetric matrices. A workaround is based on
the usual embedding

    a+ib↦[[  a  b; -b  a ]]

of the reals into complex matrices. As such, one computer the signature
of

    [[  (M)  (M); -(M)  (M) ]]

and divides by two. See in the supplementary
files <cit.>.








Certainly the resulting algorithm seems fast. Figure <ref>
shows how the time to compute a single localizer index grows at a
smaller order than the growth we saw for the Bott index. It is not
so hard to compute a single localizer index for a system with L≈1000
<cit.>. However, this is not a thorough analysis if
we are looking at a disorder averaged index study. Open boundaries
and periodic boundary conditions will have different effects for the
same system size. Also, the localizer index is, at is heart, a local
index. To see the effect of growing system size, we ought to be working
with a global index. It is possible to make the localizer index behave
globally, but this conversion is not trivial. 



Unless we are dealing with a defect (like the hole in the system used
in Figure <ref>) we are going to set
(λ_1,λ_2) to (0,0), or whatever are the coordinates
of the center of the sample according to X and Y. If we have
a disordered system and are sweeping in energy, then we may set λ_3=E_F
to a value where the localizer is, according to the computer, singular.
Technically, we should not calculate the index there, but it is often
simpler to treat 0 as positive when counting eigenvalues. This
happens so infrequently as to be a negligible effect. (This is a guess.
In fact, the author knows of no study on the accuracy of calculating
the the signature of large, sparse Hermitian matrices. Most related
work in numerical linear algebra is on estimating the number of eigenvalues
in a region <cit.>, which is of no use in the
context of numerical K-theory.)

In theory, we can avoid computing the pseudospectrum for a study of
disorder averaged index. However, we need a way to set κ,
and looking at the pseudospectrum is one way to figure values for
κ. The theory of how to set κ is the subject of continuing
research <cit.>.

Some may complain that if we need to tune the localizer index to use
it, then it is not useful. However, many measurement instruments need
tuning. In particular, when using a scanning tunneling microscope
(STM) to do spectroscopy, one must determine the sample to tip gap.
Let us focus on E_F=0 and “probing” at the center of the sample,
so λ=0. For simplicity, let us
assume there is no lattice point at (0,0). If κ is very
small the localizer index is zero because then

    L_0(κ X,κ Y,H)=[[        H κ X-iκ Y; κ X+iκ Y       -H ]]≈[[  H  0;  0 -H ]]

and the matrix on the right has spectrum that is symmetric across
0. Also the localizer index is zero if κ is large because
then 

    L_0(κ X,κ Y,H)=[[        H κ X-iκ Y; κ X+iκ Y       -H ]]≈κ[[    0 X-iY; X+iY    0 ]]

and, since X+iY is normal we again find spectrum that is symmetric
across 0.


Finding the intermediate values of κ that work can be done
in several ways. One way is to look at the joint approximate eigenvectors
<cit.> for X, Y and H that
can be extracted from the localizer, comparing the deviation of that
state in position to system size and the deviation of that state in
energy compared to the predicted bulk gap. Here we take a more visual
approach.

First we look at a slice of the pseudospectrum at two energy levels,
across a small sample, shown in Figures <ref>
and <ref>. We see that for large κ, the
image looks somewhat like STM microscopy imagery, just showing us
where the lattice points are. That is, H is almost completely ignored.
For small κ, it is the position information that is getting
ignored. What we are looking for is a clear distinction between edge
effects at E_F=0 and bulk effects at E_F=-2.4. This is all
done with the Hamiltonian used in Section <ref>,
and with L=20. Looking at Figures <ref>
and <ref> we might think 0.05≤κ≤1
can work.

We also want a stability against disorder. With the rather strong
disorder as used in Section <ref>, and just
one sample, the same slices are shown in Figures <ref>
and <ref>. At E_F=0, as shown in Figure <ref>,
we see still a clear boundary and lots of area inside the boundary
which it is safe to assume has index 0 where defined, at least
for 0.1≤κ≤1. The data in Figure <ref>
is harder to interpret without the K-theory data. What we expect,
based on the average Bott index study as in Figure <ref>,
is that the slice at E_F=-2.4 with this level of disorder should
have, approximately by area, equal parts index 0 and index 1.



Computing multivariate pseudospectrum has not caught the attention
of numerical analysts. The method the author has developed is rather
crude. It does have one optimization, which is that when a large gap
(meaning

    gap_λ(M_1,M_2,M_3)=‖(L_λ(κ X,κ Y,H))^-1‖ ^-1

or the smallest absolute value of an eigenvalue of the localizer)
is found, there is no need to compute the gap a near locations. This
is based on <cit.>. In regions where
the gap stays enough above zero (say 10^-8) where we can trust
the numerical calculations, the index cannot change. As such, well-optomized
software would first compute the pseudospectrum, ignoring regions
where the gap is large, then compute the contiguous gapped regions,
then compute the index at one point in each region. For now, we don't
have such an optomization, and compute the index everywhere where
the gap is neither too small to be meaningful or too large to be necessary.
On hopes an applied mathematician will find this of interest in the
near future.

The index data computed is as shown in the top-left panel of Figure <ref>.
The author shamelessly used image processing software to extend the
index data to contiguous regions and then used a darking-only overlay
of the index color data over the pseudospectrum. Regions that are
too dark to distinguish hue are where the local index is nearly meaningless.
It is differing index in locations with a decent gap that are of interest,
following the logic of <cit.>. This
seems like a good way to visualize the index and pseudospectrum together.

Recall that the pseudospectrum is a grey scale map on three space,
and it is regions of 3-space we should be coloring. Producing and
displaying such information is difficult to do with current algorithms.
An attempt to display the full Clifford pseudospectrum was made in
Figures <ref> and <ref>,
and one sees the large interior space, in the shape of a distorted
cube in Figure <ref> and a distorted solid
torus in Figure <ref>, where the index
is nonzero. There seem to be smaller regions that most likely have
nonzero index at higher and lower energies, but marking this by index
will be difficult to even display. 

For now, we are just trying explain what are good values of κ
for systems here of L=20. Figure <ref>
shows the index data with the pseudospectral data for κ=1
and a clean system. Figure <ref> shows the same
but now with disorder, showing “in the bulk” mainly index 1 at
E_F=0 and a blend of index 0 and index 1 at E_F=-2.4.

The fact that the left panel of Figure <ref>
shows some index 0 well in from the edge indicates a difficulty
we have in dealing with open boundaries. Changing κ does not
make this go away. What works is moving to a larger system, say L=60.
This is above where we can easily make pictures. What was advocated
in <cit.> was to keep κ constant as L
increases to create a local invariant, and to set κ=C/L for
some fixed C to create a global index. Setting κ=C/L may
well work in a gapped system. It seems that when dealing with a mobility
gap the situation can be more subtle.



§ LOCALIZER INDEX GOES GLOBAL 






We start with doing the disorder average study with L=20,40,60,80
using the localizer index, setting κ=C/L. We have an idea,
based on setting κ for L=20, that we want 0.1≤C/20≤1
and so 2≤ C≤20. Unfortunately, various values in and around
C fail to lead to convincing results, as shown in Figure <ref>.
None of these look at all like Figure <ref>. 

Looking at the plots in Figure <ref>, and
many others, it seemed the steeper plots correspond to smaller κ,
but not for all L. So perhaps it is κ→0 here
that plays the role of L→∞ when working with the
Bott index. Various calculations, for example in <cit.>,
tell us that the localizer gap often converges when L→∞.
See Lemma <ref>.

Figure <ref> shows what happens for various
values of κ, computing the disorder averaged localizer index
at a range of Fermi levels, as L increase somewhat. We emphasize
that we stick with λ_1=λ_2=0 unless we have a good
reason, like a hole in our system, so that is the one choice we make
in this section and the last. These are too computationally intensive
to extend much further, but we have soft evidence that these plots
have essentially converged at the higher values for L shown in
each instance. 



There is some theory when one has a bulk gap of a known size, as in
<cit.>. However, we have just a mobility gap.
Presumably, one could work in an appropriate C^*-algebra, find
a gapped Hamiltonian that lives in that C^*-algebra and can extend
the methods of <cit.> to that setting. For now,
we are happy to gather numerical evidence.

In Figure <ref> shows the plots for various
decreasing κ. Notice they are decreasing by powers of ten.
The plots in Figure <ref> were used to
deduce a reasonable values of L to go with each value of κ.
It is perhaps the case that open boundaries really interact with disorder
and this is what forced us into such large systems to see plots to
create Figure <ref>.



§ LOCAL NATURE OF THE LOCALIZER 


Better algorithms are needed to make the spectral localizer effective
for three-dimensional systems, and indeed would help with two-dimensional
systems. Here is one possible improvement.

Suppose one is working on a two-dimensional system with L≥300
for a specific λ. It should be possible to compute
the gap and index of a further truncation of the system to smaller
size and use those results if that gap is large enough, and only if
this is small recomputing the gap and index on the full system. This
should be possible due to the local nature of the localizer, where
drastic alterations of the Hamiltonian away from the (λ_1,λ_2)
have little effect, as we will soon see. Also, it is the smaller values
of the localizer gap that need to be computed precisely.

A glance of Figures <ref> and <ref>
shows that, at least in one instance, something drastic like cutting
out a hole in the middle of a sample has little effect on the pseudospecturm
at the edge. An analytic expression of this fact is the following
lemma. This was proving jointly with Schulz-Baldes years ago, but
in this form it never found a place in our papers. Many similar, and
more technical results, are in our papers <cit.>.

For simplicity, we state this regarding computing the spectral gap
for the localizer when λ=0, and
only for three matrices (so two physical dimensions). So let L=L_0.
Notice also the lemma is stated in terms of changes to the square
of the localizer gap.

 Suppose X, Y and H are Hermitian
matrices, with X and Y commuting, and set 

    γ=‖ L(X,Y,H)^-1‖ ^-2.

Let Z=X+iY. If H_0 is Hermitian with

    ‖ |Z|^-1(HH_0+H_0H+H_0^2)|Z|^-1‖≤ C

and

    ‖[Z,H]‖≤ D

and 

    ‖ |Z|^-1[Z,H_0]|Z|^-1‖≤ E

then

    |γ-γ_0|≤(C+E)γ+(C+E)D

where

    γ_0=‖ L(X,Y,H+H_0)^-1‖ ^-2.

If 

    γ>(C+E)D/1-(C+E)

then the index at zero for (X,Y,H+H_0) will be the same as for
(X,Y,H).



Notice Z is normal, and compute: 

    L(X,Y,H)     =[[   H Z^†;   Z  -H ]],


    L(X,Y,H)^2    =[[ H^2+|Z|^2  -[H,Z]^†;    -[H,Z] H^2+|Z|^2 ]],


    L(X,Y,H+H_0)^2    =L(X,Y,H)^2+[[ HH_0+H_0H+H_0^2      -[H_0,Z]^†;        -[H_0,Z] HH_0+H_0H+H_0^2 ]].

We deduce from (<ref>) that 

    -C|Z|^2≤ HH_0+H_0H+H_0^2≤ C|Z|^2.

From (<ref>) we deduce

    [[ -E|Z|     0;     0 -E|Z| ]]≤[[          0 -[H_0,Z]^*;   -[H_0,Z]          0 ]]≤[[ E|Z|    0;    0 E|Z| ]].

We find 


    L(X,Y,H+H_0)^2
        =[[ H^2+|Z|^2  -[H,Z]^*;    -[H,Z] H^2+|Z|^2 ]]+[[ HH_0+H_0H+H_0^2      -[H_0,Z]^*;        -[H_0,Z] HH_0+H_0H+H_0^2 ]]
       ≥[[ (1-C-E)H^2+|Z|^2         -[H,Z]^*;           -[H,Z] (1-C-E)H^2+|Z|^2 ]]+[[ (-C-E)|Z|^2           0;           0 (-C-E)|Z|^2 ]]
        =(1-C-E)B(X,Y,H)^2+[[            0 (C+E)[H,Z]^*;   (C+E)[H,Z]            0 ]]
       ≥(1-C-E)B(X,Y,H)^2-(C+E)D

and

    L(X,Y,H+H_0)^2
        =[[ H^2+|Z|^2  -[H,Z]^*;    -[H,Z] H^2+|Z|^2 ]]+[[ HH_0+H_0H+H_0^2      -[H_0,Z]^*;        -[H_0,Z] HH_0+H_0H+H_0^2 ]]
       ≤[[ (1+C+E)H^2+|Z|^2         -[H,Z]^*;           -[H,Z] (1+C+E)H^2+|Z|^2 ]]+[[ (C+E)|Z|^2          0;          0 (C+E)|Z|^2 ]]
        =(1+C+E)B(X,Y,H)^2+[[             0 (-C-E)[H,Z]^*;   (-C-E)[H,Z]             0 ]]
       ≤(1+C)B(X,Y,H)^2+(C+E)D.





§ ACKNOWLEDGMENTS


This material is based upon work supported by the National Science
Foundation under DMS 1700102. The author thanks Liang Du for pointing
out the differing use of the constants in the BHZ model. Most of the
computing was done on machines at the Center for Advance Research
Computing at the University of New Mexico. The author is most grateful
to Hermann Schulz-Baldes for allowing Lemma <ref>
to appear here.

10

Matlab_LDL
Matlab function ldl documentation.
mathworks.com/help/matlab/ref/ldl.html.
Accessed: 2019-07-24.

MUMPS_guide
MUltifrontal Massively Parallel Solver Users guide.
mumps.enseeiht.fr/doc/userguide_5.2.1.pdf.
Accessed: 2019-07-24.

supplement
Supplementary files.
math.unm.edu/%7Eloring/BottLocalizerGuide/.

bandres2016topological
Miguel A Bandres, Mikael C Rechtsman, and Mordechai Segev.
Topological photonic quasicrystals: Fractal topological spectrum and
  protected transport.
Physical Review X, 6(1):011016, 2016.

BellissardEtAl_NCG_QuantHall
J. Bellissard, A. van Elst, and H. Schulz-Baldes.
The noncommutative geometry of the quantum Hall effect.
J. Math. Phys., 35(10):5373–5451, 1994.
Topology and physics.

benzi2013decay
Michele Benzi, Paola Boito, and Nader Razouk.
Decay properties of spectral projectors with applications to
  electronic structure.
SIAM review, 55(1):3–64, 2013.

berenstein2012matrix
David Berenstein and Eric Dzienkowski.
Matrix embeddings on flat ℝ^3 and the geometry of
  membranes.
Physical Review D, 86(8):086001, 2012.

bianco2011mapping
Raffaello Bianco and Raffaele Resta.
Mapping topological order in coordinate space.
Physical Review B, 84(24):241106, 2011.

Blackadar-shape-theory
Bruce Blackadar.
Shape theory for C∗-algebras.
Math. Scand., 56(2):249–275, 1985.

ConnesHigsonAsymptoticMorphisms
Alain Connes and Nigel Higson.
Déformations, morphismes asymptotiques et K-théorie
  bivariante.
C. R. Acad. Sci. Paris Sér. I Math., 311(2):101–106, 1990.

AlessioRigolFloquetChern
Luca D'Alessio and Marcos Rigol.
Dynamical preparation of Floquet Chern insulators.
Nature communications, 6:8336, 2015.

EilersLoringContingenciesStableRelations
Søren Eilers and Terry A. Loring.
Computing contingencies for stable relations.
Internat. J. Math., 10(3):301–326, 1999.

ELP-pushBusby
Søren Eilers, Terry A. Loring, and Gert K. Pedersen.
Morphisms of extensions of C *-algebras: pushing forward the
  Busby invariant.
Adv. Math., 147(1):74–109, 1999.

ExelLoringInvariats
Ruy Exel and Terry A. Loring.
Invariants of almost commuting unitaries.
J. Funct. Anal., 95(2):364–376, 1991.

fulga2016aperiodic
Ion C Fulga, Dmitry I Pikulin, and Terry A Loring.
Aperiodic weak topological superconductors.
Physical review letters, 116(25):257002, 2016.

HastLorTheoryPractice
Matthew B. Hastings and Terry A. Loring.
Topological insulators and C *-algebras: Theory and numerical
  practice.
Ann. Physics, 326(7):1699–1759, 2011.

jiang2009numerical
Hua Jiang, Lei Wang, Qing-feng Sun, and XC Xie.
Numerical study of the topological anderson insulator in HgTe/CdTe
  quantum wells.
Physical Review B, 80(16):165316, 2009.

kisil1996mobius
Vladimir Kisil.
Möbius transformations and monogenic functional calculus.
Electron. Res. Announc. Math. Sci, 2(1):26–33, 1996.

kitaev2006anyons
Alexei Kitaev.
Anyons in an exactly solved model and beyond.
Annals of Physics, 321(1):2–111, 2006.

lin2016approximating
Lin Lin, Yousef Saad, and Chao Yang.
Approximating spectral densities of large matrices.
SIAM review, 58(1):34–65, 2016.

LoringSchuBa_even_AI
Terry Loring and Hermann Schulz-Baldes.
The spectral localizer for even index pairings.
J. Noncommut. Geom., to appear.
arXiv preprint arXiv:1802.04517.

Loring-K-thryAsymCommMatrices
Terry A. Loring.
K-theory and asymptotically commuting matrices.
Canad. J. Math., 40(1):197–216, 1988.

LoringCstarRelations
Terry A. Loring.
C^*-algebra relations.
Math. Scand., 107(1):43–72, 2010.

LoringPseudospectra
Terry A. Loring.
K-theory and pseudospectra for topological insulators.
Ann. Physics, 356:383–416, 2015.

loring2018bulk
Terry A Loring.
Bulk spectrum and k-theory for infinite-area topological quascrystal.
arXiv preprint arXiv:1811.07494, 2018.

LorHastHgTe
Terry A. Loring and Matthew B. Hastings.
Disordered topological insulators via C *-algebras.
Europhys. Lett. EPL, 92:67004, 2010.

L_S-B_finite_vol
Terry A. Loring and Hermann Schulz-Baldes.
Finite volume calculation of K-theory invariants.
New York J. Math., 23:1111–1140, 2017.

loringSchuba2019spectral
Terry A Loring and Hermann Schulz-Baldes.
Spectral flow argument localizing an odd index pairing.
Canadian Mathematical Bulletin, 62(2):373–381, 2019.

MitchellNashAmorphous
Noah P Mitchell, Lisa M Nash, Daniel Hexner, Ari M Turner, and William TM
  Irvine.
Amorphous topological insulators constructed from random point sets.
Nature Physics, 14(4):380, 2018.

ProdanTopInsNCG
Emil Prodan.
Disordered topological insulators: a non-commutative geometry
  perspective.
J. Phys. A, 44(11):113001, 50, 2011.

sykora2016fuzzy
Andreas Sykora.
The fuzzy space construction kit.
arXiv preprint arXiv:1610.01504, 2016.

toniolo2018time
Daniele Toniolo.
Time-dependent topological systems: A study of the bott index.
Physical Review B, 98(23):235425, 2018.

V_S_S-B_HalfSignature
Jonas; Viesca, Edgar Lozano; Schober and Hermann Schulz-Baldes.
Chern numbers as half-signature of the spectral localizer.
J. Math. Phys., to appear.



