Open boundary condition in hubbard model In case of open boundary condition, crystal momentum is not a good quantum number anymore. I mean we can construct an $N×N$ matrix in position basis and diagonalize it to obtain the eigenvalues, but how do we diagonalize it analytically? That is, how do we obtain band structure in such cases?
 A: Well, first of all the Hamiltonian of the Hubbard model has both an hopping term and an interacting term of the form
$$
H = -t \sum_{ij,\, \sigma} ( c^{\dagger}_{i\sigma} c_{j\sigma} + \mathrm{h.c.} ) + U \sum_i n_{i\uparrow} n_{i\downarrow},
$$
where $n_{i\sigma}$ is the number operator associated to the fermionic operators $c^{(\dagger)}_{i\sigma}$.
I think you are interested in the simple case $U=0$, where the Hamiltonian is quadratic, the particles are non-interacting and the whole problem boils down to compute the energy states available for a single particle.
In this case, if you have periodic boundary conditions (PBC) on the lattice, you can switch to the momentum representation (Fourier transform the fermionic operators) and get the band structure.
If you have open boundary conditions (OBC) on the lattice, momentum is no longer a good quantum number and you have to diagonalize an $N\times N$ matrix ($N$ being the number of sites) written in the Fock basis of a single localized fermion (as you said).
As far as I know there is no simple solution that you can write on a piece of paper for the eigenstates when $N$ is very large, but a numerical diagonalization is relatively fast and gives you all the energy states quite effectively.
When $U\neq 0$ the story is completely different. In this case in principle you have to write a $4^N \times 4^N$ matrix (both for PBC and OBC), where $4$ is the local number of fermionic configurations (empty site, spin up particle, spin down particle, double occupation) and $4^N$ is the dimension of the whole Hilbert space of many body Fock states.
This calculation is much more difficult when $N$ is large, because the matrix size scales exponentially with the lattice size. There are many strategies to simplify the diagonalization that if you want I can discuss, but the exponential scaling remains a big issue.
A: Most likely you are thinking about Hubbard model with $U=0$. So we really just have a hopping model:
$$
H=-\sum_{j=1}^{N} (c_j^\dagger c_{j+1}+\text{h.c.})
$$
The spin indices are suppressed because without interaction they are completely decoupled. Here I assume it's 1-dimensional system.
With open boundary conditions, the matrix in the position space we want to diagonalize becomes
$$
h=-\begin{pmatrix}
    0 & 1 & & & 0 \\
    1 & 0 & 1 & & \\
    & \ddots & \ddots & \ddots & \\
    & & 1 & 0 & 1 \\
    0 & & & 1 & 0
    \end{pmatrix}
$$
Let us look for the eigenvectors. Denote the eigenvector by $v=(v_1,v_2,\cdots, v_N)$, and the energy $E$. Then we need to solve
$$
v_{j-1}+Ev_j+v_{j+1}=0, j=1,2,\cdots,N
$$
Here I define $v_0=v_{N+1}=0$.
It is easy to see that the solutions are $v_j=\sin kj$, where $k$ satisfies $\sin (N+1)k=0$, so $k=\frac{\pi}{N+1}n$ for $n=1,2,\cdots,N$, with energy $E=-2\cos k$.
