Tight binding in the limit of large system size Suppose that one has a continuous Hamiltonian with spin-orbit interaction, for example
$H=-\dfrac{\mathbf{p}^2}{2m} +\kappa({\boldsymbol\sigma}\times\mathbf{p}) + U(x)$
and want to approximate this Hamiltonian with a tight-binding discrete model.
This can be done with the method of finite differences (e.g., Datta, Quantum Transport), using the substitution
$(\partial^2_x \psi)_{x=x_n} \rightarrow \frac1{a^2}[\psi(x_{n+1}-2\psi(x_n)+\psi(x_{n-1})]$
$(\partial_x \psi)_{x=x_n} \rightarrow \frac1{2a}[\psi(x_{n+1}-\psi(x_{n-1})]$
Suppose now that one has the opposite problem, i.e., one has a tight-binding Hamiltonian and want to obtain informations on the system in the continuous limit.
Just as an example, consider the tight-binding Hamiltonian
$H=\sum_{i=1}^N\sum_{ss'} \delta_{ss'} [u(r_i)-\mu] a^\dagger_{is} a_{is} + \delta_{ss'} t a^\dagger_{is} a_{i+1,s} +  \imath\alpha\sigma^y_{ss'} a^\dagger_{is} a_{i+1,s'} + \text{h.c.}$
where $s$ and $t$ are the spin indexes.
This Hamiltonian is in the general case not analytically solvable, since the parameters and $u(r_i)$ can depend on the lattice sites (e.g., the system can contain one or more impurities).
However, for finite $N$ one can always, in principle, diagonalize the Hamiltonian numerically and calculate some relevant physical quantities, for example the gap.
I am explicitly interested in cases where the tight-binding model can be solved only numerically.
Is there a general method to take the continuous limit $N\rightarrow\infty$ and $Na=\text{constant}$ of a general discrete tight-binding model?
To avoid confusion, I stress that I am not interested to the case where the system is infinite, but rather on the limit case where the system is continuous, i.e., not discrete, but of finite size $L=Na$.
One idea is to calculate numerically the spectra and the relevant physical quantities for increasing $N$ and with a decreasing lattice parameter $a\rightarrow 0$ with the constraint $Na=\text{constant}$ and to deduce the asymptotic behavior.
Since in the continuous limit $t=\hbar^2/(2m a^2)$ and $\alpha=\kappa/(2a)$ ($\kappa$ is the SO parameter), this is equivalent to take the limit $N\rightarrow\infty$ and taking $t\propto N^2$ and $\alpha\propto N$.
Is this approach correct?
Is it correct to assume the limit $N\rightarrow\infty$ with $t\propto N^2$ and $\alpha=\propto N$  at each step of the calculation?
Is there a reference or a book that discuss the discrete to continuous limit?
Many thanks in advance.
 A: To obtain the energy gap in the thermodynamic limit, one should take $N\rightarrow\infty, V\rightarrow\infty$ where $N$ is the number of atoms and $V$ the system size, but hold $N/V$ (i.e. density) fixed. In your case, it means simply taking $N$ to $\infty$ is enough and keep all $t, u, \alpha$ fixed for now. This tight-binding model can be solved exactly. First of all, observe that the first two terms are identity if the spin part is concerned, so one is free to choose any spin quantization axe. Let us for simplicity just rotate the spin axes such that the last term is proportional to $\sigma_z$. Assuming periodic boundary condition $a_{j+N}=a_j$, you can go to the momentum space $a_j=\frac{1}{\sqrt{N}}\sum_k a_k e^{ikj}$, the Hamiltonian becomes
$H=\sum_k (2t\cos k+u)a_{ks}^\dagger a_{ks}+2\alpha\cos k(a_{k\uparrow}^\dagger a_{k\uparrow}-a_{k\downarrow}^\dagger a_{k\downarrow})$
Basically spin up and down electrons are decoupled. Spin up/down electrons have a dispersion
$E_{\uparrow/\downarrow}(k)=(2t\pm\alpha)\cos k+u$. So there are two bands, and now you can just fill up the bands with the electrons.
Notice that $u$ is (minus) the chemical potential, which tells you basically energy levels should be filled up to $-u$. Notice that the bands for spin up and down electrons have almost the same "shape", so you can easily convince your self that as long as $-u$ crosses the dispersion at all (meaning that the bands are neither empty nor fully occupied), you will end up with a metal which has no gap in the thermodynamic limit.
Just one more comment:I'm guessing you are trying to model 1d electrons with spin-orbit coupling. Usually, in systems like semiconductor nanowires, the spin-orbit coupling (e.g. Rashba) looks like $ia_i^\dagger \sigma_y a_{i+1}+\text{h.c.}$ (notice the $i$ in front), not the one you wrote.
