Wannier sectors

[bfocassio]

Can you explain how does the pw2z2pack code perform the separation of Wannier eigenstates into two symmetry sectors or give me any reference to understand how to perform such separation? (The physicis behind it) How one can separate the Wannier eigenstates into two sectors for computing the mirror Chern number or spin chern number?

Is it possible to do such thing using Hamiltonians in z2pack?

[greschd]

Sure! The physical background is that you have a symmetry that maps k into itself - for some subset of the Brillouin zone, for example a mirror plane. Because the Hamiltonian and symmetry commute, you can simultaneously diagonalise them - in other words, the Hamiltonian is block-diagonal w.r.t. the symmetry eigenspaces.

To compute the topological invariants corresponding to each symmetry eigenvalue, you just restrict the calculation to the symmetry eigenspace.

You can do this, as long as you know the symmetry representation in the basis you're using. It's all pretty nicely described in this semester thesis: semester_thesis_tmetger.pdf
The code from that thesis is in the dev/symmetry branch - unfortunately I haven't gotten around to fully integrating / releasing it.

[bfocassio]

Thank you very much for the answer. Reading this paper really helped.

I'm trying this with a ''simple'' example of calculating the spin Chern number for the Kane-Mele model.

As found in the examples folder, we can define the following:

import logging

import z2pack
import numpy as np
from numpy import cos, sin, kron, sqrt
import matplotlib.pyplot as plt

logging.getLogger('z2pack').setLevel(logging.WARNING)
import scipy.linalg as la

# defining pauli matrices
identity = np.identity(2, dtype=complex)
pauli_x = np.array([[0, 1], [1, 0]], dtype=complex)
pauli_y = np.array([[0, -1j], [1j, 0]], dtype=complex)
pauli_z = np.array([[1, 0], [0, -1]], dtype=complex)

def get_kane_mele_hamiltonian(t, lambda_v, lambda_R, lambda_SO):
    def inner(k):
        k = np.array(k) * 2 * np.pi
        kx, ky = k
        # change to reduced coordinates
        x = (kx - ky) / 2
        y = (kx + ky) / 2
        return (
            t * (1 + 2 * cos(x) * cos(y)) * kron(pauli_x, identity) +
            lambda_v * kron(pauli_z, identity) + lambda_R *
            (1 - cos(x) * cos(y)) * kron(pauli_y, pauli_x) +
            -sqrt(3) * lambda_R * sin(x) * sin(y) * kron(pauli_y, pauli_y) +
            2 * t * cos(x) * sin(y) * kron(pauli_y, identity) + lambda_SO *
            (2 * sin(2 * x) - 4 * sin(x) * cos(y)) * kron(pauli_z, pauli_z) +
            lambda_R * cos(x) * sin(y) * kron(pauli_x, pauli_x) +
            -sqrt(3) * lambda_R * sin(x) * cos(y) * kron(pauli_x, pauli_y)
        )

    return inner

Now the symmetry part should be something as:

def symmetry():
    return la.block_diag(*([pauli_z] * 2))

Then onde could simply do:

system = z2pack.hm.System(
    get_kane_mele_hamiltonian(t=1, lambda_v=0.1, lambda_R=0.05, lambda_SO=0.06),
    dim=2,
    check_periodic=True,
    symm=symmetry())

res = z2pack.surface.run(system=system,
                         surface=lambda s, t: [s, t],
                         use_symm=True)

C_plus = z2pack.invariant.chern(res.symm_project(1))
C_minus = z2pack.invariant.chern(res.symm_project(-1))
C_plus,C_minus,C_plus-C_minus

This gives me the following :

(1.4936833116537378, 1.0, 0.4936833116537378)

which should be integers

Am I missing something in this representation of the symmetry ?

[greschd]

I think the problem here is that time-reversal maps k onto -k. So on any surface you choose, it will not map a given k back into itself.

Note that the method actually relies on this: at each k, the Hamiltonian needs to commute with the symmetry to become block-diagonal. I think the same is true for pw2z2pack - but I'd have to check the source again, since that code was written by Gabriel Autès.

In principle it is possible to separate the Hilbert space according to symmetry eigenvalues even if they are not local in k, but that isn't implemented here (it would be a "global" operation, if you will). But again, time-reversal is a special case: We put a very short explanation of that in Appendix F of the Z2Pack paper.

[bfocassio]

Thank you very much

[bfocassio]

Sorry to reopen this.

I'm very interested in the symmetry projection feature and trying to do a similar implementation on other codes that use Quantum Espresso projections on atomic wavefunctions.

The remaining question is: Given a calculation with SOC in QE and the projection in atomic wavefunctions (as given by projwfc.x), how can I write the mirror operator ? or any other rotation operator.

My attempt is: Since projwfc gives me the wfc on the |j,m> basis, I should write the angular momentum operators J_x, J_y, and J_z for a given j and m.
Then I would obtain the rotation matrix by D = exp(-i*theta*J).
The final operator is constructed by arranging the D matrices for each j at the diagonal given the same order as in projwfc output.

For example, consider a system with mirror symmetry in the y-axis, the M_y operator (mirror y) would be given by using J_y and theta = pi for building the D matrices. For a system with wfc order as j=0.5, j=1.5 this would be M_y = diag(D_0.5,D_1.5).
Is this correct ?

Is there any special consideration regarding the |j,m> basis in QE ?

I would appreciate any help.

[greschd]

I'm not entirely sure if it's what you're trying to do, but the symmetry-representation code can generate representation matrices from the symmetry + orbital information.

The only special consideration I can think of is |j,m> not being an orthogonal basis, but I think you're already aware of that (#79).

[bfocassio]

Thank you very much, the code is truly helpful. I managed to construct the matrices.

[greschd]

Fantastic! Out of curiosity: did you end up using the https://github.com/Z2PackDev/Z2Pack/tree/dev/symmetry branch? If so, how did it work for you?

[bfocassio]

Sorry for the late response

I'm testing the symmetry feature with M. Ezawa's TCI model with C4 symmetry.

This TCI phase is protected by mirror-z symmetry which the author writes as

Setting up the Pauli matrices as:

eye = np.eye(2,dtype=complex)
sx = np.array([[0,1],[1,0]],dtype=complex)
sy = np.array([[0,-1j],[1j,0]],dtype=complex)
sz = np.array([[1,0],[0,-1]],dtype=complex)

To implement the model I use:

def do_Hk(k,m0,m1,l1,l2,Ez):

    from scipy import sparse

    k = np.array(k) * np.pi * 2
    kx,ky,kz = k[0],k[1],k[2]
    
    Ax = l1 * np.sin(kx) + l2 * np.sin(kx) * np.cos(ky)
    Ay = l1 * np.sin(ky) + l2 * np.sin(ky) * np.cos(kx)
    
    m = m0 + m1 * (np.cos(kx) + np.cos(ky))
    
    Hm = sparse.kron(m*eye,sx).toarray()
    
    Hso = sparse.kron(Ax*sx,sy).toarray() + sparse.kron(Ay*sy,sy).toarray()
    
    Hz = Ez * sparse.kron(eye,sz).toarray()

    Hk = Hso + Hm + Hz

    return Hk

The symmetry is then Mz = -1j*sparse.kron(sz,sx).toarray()

Then I wish to compute the mirror Chern number wich is given by

Therefore, we should project the result into the +i and -i eigenspaces and compute the individual Chern numbers.

This could be done by:

m0=-0.6
m1=0.4
lambda_1=1.0
lambda_2 = -0.5
Ez=0.0

system = z2pack.hm.System(lambda k: do_Hk(k,m0=m0,m1=m1,l1=lambda_1,l2=lambda_2,Ez=Ez), dim=3, hermitian_tol=None,symm=Mz)
result = z2pack.surface.run(system=system, surface=lambda s, t: [s,t,0.], min_neighbour_dist=1e-5,use_symm=True)

To compute the Chern numbers I use:

projected_results = [result.symm_project(i) for i in [1.j,-1.j]]
total_chern = z2pack.invariant.chern(result)
individual_chern = [z2pack.invariant.chern(r) for r in projected_results]

In the above I'm using some parameters given by the author in Fig. 3(d3) which he suggests that yields C_M = 2.

First I noticed that in the symmetry features, eigenvalues are computed using np.eigh which enforces real eigenvalues

Even after changing the diagonalization algorithm to eig and rounding the result (inside hm.py), I could not get the correct result.

Can you help me ?

[greschd]

Hi @bfocassio, thanks for the report! I'll have to dig into this a bit to see what's going on. Just to manage expectations: It'll probably take a week or two until I get around to having a closer look at this.