Momentum creation operator

I am looking at the Bose-Hubbard model. Specifically the solution in the limit of weak interactions $U = 0$. I understand the proof to show the Hamiltonian is diagonal in momentum space, however I would like clarify the origin of the transformation

$$\hat{b}_{i}^{\dagger} = \frac{1}{\sqrt{M}}\sum_{\mathbf{k}}\hat{c}_{\mathbf{k}}^{\dagger}e^{-i\mathbf{k}\cdot\mathbf{r}_{i}}$$

where $\hat{b}_{i}^{\dagger}$ is the creation operator for site $i$. $\mathbf{k}$ is the momentum vector, $\mathbf{r}_{i}$ is the position vector of the $i$th lattice site. $\hat{c}_{\mathbf{k}}^{\dagger}$ is the momentum creation operator. $M$ is the number of lattice sites.

Several questions. Firstly I understand the sum $\sum_{\mathbf{k}}$ is a sum over momentum space which is discrete because the position space is discrete, why does a discretised position space imply a discretised momentum space?

Secondly I understand this equation is a sort of fourier transform between position and momentum space operators. But it can also be determined by changing the basis of the creation operator, could somebody show me how this is explicitly done?

I guess I am struggling to understand this dual space picture (between position and momentum).

Lastly $$\sum_{i = 1}^{i = M}e^{i(\mathbf{k}^{'} - \mathbf{k})\cdot\mathbf{r}_{i}} = M\delta_{\mathbf{k},\mathbf{k^{'}}}$$ where $\mathbf{k}^{'}$ is another momentum vector. Is the orthonormality of this sum because $e^{i\mathbf{k}\cdot\mathbf{r}_{i}}$ is eigenfunction of the momentum operator and so it is orthogonal to any other eigenfunction $e^{i\mathbf{k^{'}}\cdot\mathbf{r}_{i}}$. Does this plane wave form some sort of of basis for both position and momentum space or is that wrong and it is just representative of the overlap between the spaces?

A discrete position space does not imply a discrete momentum space. A discrete positions space implies a bounded momentum space, and a bounded position space implies a discrete momentum space. To see this, you should consider a number of points, $N$, arranged on a ring. We can label sites by an index $n$, with the property that $n=N$ is the same as $n=0$. If we have some function (or operator) defined on these points, say $a_n$, we want to write out the function as a Fourier series:
$$a_n=\sum_k a_ke^{ikn}$$
What properties do the $k$ need to have? First, we know $a_n$ is periodic, so $a_N=a_0$. This implies that $\sum_k a_ke^{ikN}=\sum_k a_ke^{ik0}$, which only holds if $kN$ is a multiple of $2\pi$. Thus, we must have $k=\frac{2\pi m}{N}$ for some integer $m$. Thus, we've seen that the fact that $N$ is finite implies $k$ is discrete (and if we let $N\rightarrow \infty$, $k$ would become continuous).
Now, we use the fact that $a_n$ is discrete. Because n is an integer, we have that $e^{ikn}=e^{i(k+2\pi)n}$ for any $k$. That means we only care about $k$ modulo $2\pi$. In other words, we can choose all our $k$'s to fall between $\pi$ and $-\pi$.