Numerical Berry curvature for bosons

If $|n(mathbf{k})rangle := alpha^dagger_mathbf{k}|0rangle$, then the Berry curvature is undefined: that's the dilemma that was highlighted in the question. The states $|n(mathbf{k})rangle$ typically used in Berry-curvature calculations cannot be the states $alpha^dagger_mathbf{k}|0rangle$ that were defined in the question.

The states $|n(mathbf{k})rangle$ typically used in Berry-curvature calculations (illustrated below) are eigenvectors of a matrix $H(mathbf{k})$ that is distinct from the operator that was denoted $E(mathbf{k})alpha_mathbf{k}^daggeralpha_mathbf{k}$ in the question, but they are related: the coefficient $E(mathbf{k})$ is an eigenvalue of the matrix $H(mathbf{k})$ as well as being an eigenvalue of the operator $E(mathbf{k})alpha_mathbf{k}^daggeralpha_mathbf{k}$. However, the eigenvectors $|n(mathbf{k})rangle$ of $H(mathbf{k})$ live in a different Hilbert space from the eigenstates of $E(mathbf{k})alpha_mathbf{k}^daggeralpha_mathbf{k}$.

To illustrate this, consider equation (6.45) in [1]: $$ newcommand{da}{a^dagger} newcommand{db}{b^dagger} newcommand{la}{langle} newcommand{ra}{rangle} newcommand{pl}{partial} H = sum_mathbf{k} Big(da(mathbf{k}) db(mathbf{k})Big) H(mathbf{k})left(begin{matrix} a(mathbf{k}) b(mathbf{k})end{matrix}right) tag{1} $$ where $H(mathbf{k})$ is a $2times 2$ matrix. The Hamiltonian (1) can be derived from a special case of the Hamiltonian shown in the question, using appropriately-defined Fourier transforms.$^{[2]}$ This model has two bands.$^{[3]}$ In this case, the states $alpha^dagger_mathbf{k}|0ra$ described in the question are linear combinations of $da(mathbf{k})|0ra$ and $db(mathbf{k})|0ra$, where $|0ra$ is the ground state of (1), and the coefficients in this linear combination are the components of one the eigenvectors of the matrix $H(mathbf{k})$. Different eigenvectors of $H(mathbf{k})$ correspond to different bands. The states $alpha_mathbf{k}^dagger|0ra$ are mutually orthogonal, as observed in the question, so we cannot define a Berry curvature for them. In contrast, the states $|n(mathbf{k})ra$ typically used in Berry-curvature calculations do not belong to the original Hilbert space at all. They belong to a different Hilbert space that is two-dimensional (in this two-band example), namely the Hilbert space in which the matrix $H(mathbf{k})$ is defined. If we take the state $|n(mathbf{k})ra$ to be one of the eigenstates of this $2times 2$ matrix, parameterized by $k$, then the states in this one-parameter family are typically not orthogonal to each other, so they can have a well-defined Berry curvature.

If this is the intent of the question, then the Berry curvature can be computed numerically using the approach in cond-mat/0503172, as reviewed in the question.

[1] http://www-personal.umich.edu/~sunkai/teaching/Fall_2014/Chapter6.pdf

[2] For the model considered in [1], a Bogoliubov transform is not needed, but it can be generalized so that a Bogoliubov transform is needed.

[3] The model can be generalized to an arbitrary number of bands, in which case $H(mathbf{k})$ is generalized to a matrix of arbitrary size.

There are multiple things to address here. First you would have reached the same conclusion whether the particles are bosons or fermions. In your reasoning as you defined $|n(textbf{k})rangle = alpha^dagger_textbf{k} |0rangle$, $langle n(textbf{k})|n(textbf{k} + delta textbf{k})rangle = 0$, whether $[alpha_textbf{k}, alpha^dagger_{textbf{k}^prime}] = delta_{textbf{k}textbf{k}^prime}$ or ${alpha_textbf{k}, alpha^dagger_{textbf{k}^prime}}=delta_{textbf{k}textbf{k}^prime}$. However I do think it's an interesting question to ask what's the meaning of the Berry curvature defined this way for a system of bosons. That's a different question though.

Before getting to the main confusion, I also wanna stress that you ran into a numerical problem. Mathematically you haven't shown that the $U(1)$ link is zero. You've shown that the numerator is zero, but you also have the denominator to be zero. So strictly speaking, one would need to work harder. However of course I understand, if you are running some sort of code this is bad.

Chiral Anomaly's answer address the confusion about the different Hilbert spaces. Here I'll try to be more explicit. The paper you cited is written in the context of condensed matter systems, that is periodic systems with some discreet transnational symmetries. For such periodic systems one always make use of Bloch's theorem and write the wavefunction as, $$psi^n_textbf{k}(textbf{r}) = e^{itextbf{k} cdottextbf{r}} u^n_textbf{k}(textbf{r}),$$ where $u^n(textbf{r})$ is a function that doesn't transform under the transnational symmetry of the system. It can be easily verified that $$langle psi^n_textbf{k}(textbf{r})|psi^{n^prime}_{textbf{k}^prime}(textbf{r})rangle = delta_{textbf{k}textbf{k}^prime}delta_{nn^prime}.$$ Now in your question we have $$alpha^{n dagger}_textbf{k}|0rangle = |psi^n_textbf{k} rangle,$$ and with our choice of normalization above, this is consistent with your commutation relationships for $alpha^n_{textbf{k}}$, and $alpha^{ndagger}_{textbf{k}}$, because remember that in general. $$[alpha^n_textbf{k}, alpha^{ndagger}_{textbf{k}^prime}] = langle psi^n_textbf{k}(textbf{r})|psi^{n^prime}_{textbf{k}^prime}(textbf{r})rangle.$$ The wavefunction used in the paper you referred to, $|n(textbf{k})rangle$, is not the same. Rather with these periodic systems usually it's convenient to restrict oneself to only one unit cell. Here $|n(textbf{k})rangle$ would be defined within one unit cell of the system such that, $$langle r|n(textbf{k})rangle equiv u_textbf{k}^n(textbf{r}).$$
What I mean by within a unit cell of the system is that when we see expressions like $langle n(textbf{k}) | n(textbf{k}^prime)rangle $, the integration is only taken within one unit cell. Notice again that $| n(textbf{k})rangle$ is not what you have as $alpha^n_{textbf{k}}|0rangle $.

Now the all important point is that even though $langle psi^n_textbf{k}(textbf{r})|psi^{n}_{textbf{k} + delta textbf{k}}(textbf{r})rangle = 0$, and hence the confusion, in general $langle n_textbf{k}(textbf{r})|n_{textbf{k} + delta textbf{k}}(textbf{r})rangle neq 0$. I'm not fully aware how your code works for the calculation, but basically you need to extract the periodic part from your $alpha^dagger_{textbf{k}}|0rangle$, and then put that into the formula for the $U(1)$ link.

You might ask why use $|n(textbf{k})rangle$ in the definition of the $U(1)$ link and not $|psi_{textbf{k}}rangle$. We already mentioned the awkward zero over zero situation. But also note that as one move through any path in the Brillouin zone, the Berry phase gained by $|psi_{textbf{k}}rangle$ is the same as the berry phase gained by $|n(textbf{k})rangle$, and so it's usefull to define the $U(1)$ link in terms of $|n(text{k})rangle$.