How to numerically solve real-space 1D time-independent Schrodinger equation using 2D momentum-space Hamiltonian?

Question

Consider the usual simple 2-level graphene Hamiltonian with mass in momentum-space where:

$$ H(k,V)=-t \sum_{\delta} [\cos(k\cdot\delta)\sigma_x-\sin(k\cdot\delta)\sigma_y+V\sigma_z], $$ where $t$ is the hopping, $\sigma_i$ are Pauli matrices, $\delta$ lattice vectors and $V$ an on-site potential ($V=0$ for graphene, where the energy dispersion is degenerate/gapless at a Dirac point).

The parameter space above is 2D momentum/k-space, and the eigenvalue equation can be easily solved as I show in my answer here: How to create a 3D band structure from DFT band structure calculation?

Now, I want to numerically solve its real-space analog. I know k and real space are related to each other via Fourier transforms (https://quantummechanics.ucsd.edu/ph130a/130_notes/node82.html). However, I want to plot the wavefunction amplitudes of this system along a 1D real space line where the on-site potential $V(x)$ varies. So, I want to plot something like figure c) here: https://arxiv.org/pdf/2005.06096.pdf (reproduced below):

Is there a simple way to do this? I know I am being naive here, but would it be as simple as just taking the Fourier transforms of eigenvectors from the k-space method? This comes with the issue of choosing k-points corresponding to the real-space 1D line. I know that in theory, people usually change the Hamiltonian between spaces by inverting the creation/annihilation operators' spaces, etc. But is there something more straightforward numerically (something as silly as Fourier transforming the above k-space Hamiltonian)? There's also just replacing $k_x$ with $~\partial_x$ (as in section 3.5.7 in https://arxiv.org/abs/1310.0255), but I am just wondering whether there is something more straightforward due to numerics.

Or should I come up with some other way? I do not wish to use DFT and ab initio techniques, just simple numerical software like Python or MATLAB.

Answer 1

As you mention, your graphene hamiltonian operates on sites, i.e. it is a lattice model. Its solution therefore does not provide a wave function in real space, it provides wave function amplitudes on the graphene lattice.

That said, given the amplitudes on the graphene lattice, you can create an approximate real-space wave function from an eigenstate of the lattice hamiltonian through linear combination of atomic orbitals: Place a carbon $p_z$ orbital on every site and simply sum them up with the corresponding amplitudes (which can easily be done analytically).

For armchair and zigzag graphene nanoribbons without onsite potential the solutions of this nearest-neighbor single-orbital tight-binding hamiltonian can be obtained analytically. For playing around with different graphene nanoribbon shapes and on-site potentials I've found the pythtb code very useful (e.g. I used it to create Figure 2.5 in my Ph.D. thesis). You can use that to get the amplitudes on the lattice from the eigenstates of the hamiltonian and then do your postprocessing to get the wave function in real space.