confusion in discrete transform to solve kronig penney matrix equation in fourier space

Question

I have a periodic potential $$V(x) =\sum_{K}e^{iKx}V_{K} =\sum_{n}e^{\iota2\pi nx/a}V_{n} $$ where $K =\frac{2\pi n}a$ is the reciprocal lattice vector and $a$ is the lattice constant and $n =\pm 0,\pm 1,\pm 2,\pm 3... $so on. I would like to find the fourier coefficients $V_{K}=V_{n}$ corresponding to a particular $K$ or $n$. Suppose I have a vector for $V(x)$ having 10000 points for $$x = 0,0.01a,0.02a,...a,1.01a,....2a....99.99a$$ such that the size of my lattice is $100a$. Now $n$ will also go from $-50$ to $+49$. Thus I have defined the potential for 10000 points on a 1D lattice of 100 atoms. FFT on this vector gives 10000 Fourier coefficients. By theory of discrete fourier transform (http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf), the $K$ values corresponding to these fourier coefficients are $\frac{2\pi n}{(NX)}$ where $N$ is the no. of readings = 10000 each separated by spacing $X=0.01a$ with $n=0,1,2,3,...9999$.

But I started with $K$ that had the form $K =\frac{2\pi n}a$. What am I missing ? How do I correctly find the fourier coefficients $V_{K}$ numerically using DFT (Fast Fourier transform method in Matlab can be used) ?

For reference on the Kronig Penney fourier space matrix equation, look here http://www.physics.buffalo.edu/phy410-505/topic5/index.html

Answer 1

First, a somewhat minor point is that $x = 0,0.01a,0.02a,...a,1.01a,....2a....100a$ actually gives a list of 10001 points, not 10000 points. I will assume that you actually meant to say $x = 0,0.01a,...a,1.01a,....2a....99.99a$.

Second, you say that $$V(x)=\sum_{K}e^{iKx}V_{K}$$ where $K =\frac{2\pi n}a$ and $n=0,1,2,3$, but this gives a non-Hermitian Hamiltonian, so I will instead assume that you actually meant to say $n=0,\pm1,\pm2,\pm3$ and where $V_{-K}=V_K^*$.

So the easiest way to interpret the DFT is as follows. Note that the $k$th entry (where $k=1,2,...,10000$) of the vector $\mathbf{V}$ can be written as $$\mathbf{V}[k]=\sum_{n=-3}^3V_ne^{2\pi i 100 n (k-1)/10000}$$ where I have specifically not canceled the factor $100/10000$ in order to make the frequency domain form of the coefficients clear. From this, it's clear that each coefficient $V_n$ can be determined from the DFT as $$V_n=\sqrt{\frac{1}{10000}}\text{DFT}(\mathbf{V})[100n+1\text{ mod }10000]$$ where $\text{DFT}(\mathbf{V})$ is the unitary FFT of $\mathbf{V}$.

As a numerical example in Mathematica:

V0 = 10;
V1 = 2 + I;
V2 = 3 + 2 I;
V3 = 4;
Vm1 = 2 - I;
Vm2 = 3 - 2 I;
Vm3 = 4;
V = Chop@Table[
    Exp[I 2 \[Pi] 1/a x] V1 + Exp[I 2 \[Pi] 2/a x] V2 + 
     Exp[I 2 \[Pi] 3/a x] V3 + Exp[-I 2 \[Pi] 1/a x] Vm1 + 
     Exp[-I 2 \[Pi] 2/a x] Vm2 + Exp[-I 2 \[Pi] 3/a x] Vm3 + V0, {x, 
     0, 99.99 a, 0.01 a}];
DFTV = InverseFourier[V];

Here is the first unit cell:

ListLinePlot[V[[1 ;; 100]]]

Here is the DFT:

ListLinePlot[Abs[DFTV], PlotRange -> All, Frame -> True, 
Axes -> False]

And here are the coefficients returned from the DFT:

Chop@Table[1/Sqrt[10000] DFTV[[Mod[100 n + 1, 10000]]], {n, -3, 3}]

(*Out: {4., 3. - 2. I, 2. - 1. I, 10., 2. + 1. I, 3. + 2. I, 4.} *)

You can verify yourself that all the other DFT coefficients are zero up to machine precision. Minor note: I used InverseFourier rather than Fourier because Mathematica's definition of Fourier and InverseFourier are reversed compared to how people commonly define Fourier transforms.