A comprehensive study of the velocity, momentum and position matrix elements for Bloch states: application to a local orbital basis

We are very grateful to the referee for the thorough and extensive review work. Here we address all his/her comments in the following. We hope that the concomitant modifications we have introduced in our manuscript have turned it into a more clear and rigorous work, and also more accesible to read.

We also make a resubmission including changes motivated by our referee's suggestions, as well as other minor changes.

1. I don't really understand why the manuscript puts so much emphasis on the result of Gu et al., Eq. (1). I was unfamiliar with this result, and had a hard time understanding at first what it was about.

The work by Gu et al. deals with the commutator formula for velocity (or momentum, as they consider) following a finite volume normalization (PBCs) and a representation using coordinate space. The main conclusion there is that momentum and dipole matrix elements (multiplied by the frequency) should never be interchanged when dealing with Bloch eigenstates in a finite volume. We consider that such reference treats excelently an important part needed to understand the formalism of momentum and position operators among the vast literature in this topic. Therefore, we felt it had to be sufficiently covered in our work.

1.(a) The sentence below Eq. (1) suggested to my mind a finite crystallite ofvolume $v$ with physical surfaces, so I had trouble understanding why the eigenstates were not ``decaying at the boundaries" (vacuum tails). The authors should explain at the outset that the volume in question is just a mathematical construction of an $N \times N \times N$ supercell inscribed in the bulk. But just below Eq. (1) and later in Sec. II.D, it is clear that this expression is not true for general wave functions, or even for arbitrary functions of Bloch form; it is only true for Bloch states (eigenstates of $\hat{H}$). Also, it seems only to be true for $\hat{H}' = 0$ (absence of nonlocal terms in the potential). The authors also point out that it involves an awkward set of surface integrals. I wonder if this entire discussion could be deemphasized, or moved to an Appendix.

As pointed out by the referee, Eq. (1) is derived assuming periodic boundary conditions for eigenstates, which implies an integration domain of volume $v$ when dealing with scalar products. We have rewritten the comment below Eq. (1) in order to clarify this detail. Regarding the second point, we feel that is appropiate to mantain the discussion in Sec. II.D in the main text as we think it greatly contributes to the goal of Sec. II. This is, explaining in a comprehensive manner the different approaches to deal with the relations between the velocity, momentum and position operators.

1.(b) I have a few comments about Eqs. (11-12) and the surrounding discussion. First, isn't it true that $\langle n\boldsymbol{k}|\hat{\boldsymbol{v}}|n'\boldsymbol{k}'\rangle_v$ in Eq. (1) contains a $\langle n\boldsymbol{k}|\hat{\boldsymbol{v}}|n'\boldsymbol{k}\rangle_v$? Why don't I see this reflected in Eqs. (10-12) Shouldn't we be writing expressions for quantities like $\langle n\boldsymbol{k}|\hat{\boldsymbol{v}}|n'\boldsymbol{k}\rangle_v$ instead? Are the two terms in Eq. (11) individually nonzero in Eq. (11), but they cancel for $\boldsymbol{k} \neq \boldsymbol{k}'$? I find it confusing.

Eqs. (11) and (12) are derived by directly projecting in the coordinate space of volume $v$. With this procedure, it is not possible to factorize out a kronecker delta and directly see the diagonal behaviour. As said in the question, both terms can be individually nonzero but they cancel for transitions between states with different crystal vectors. This is explicitely shown in Ref. [4] (Gu et al.) with an example. On the other hand, this diagonal behaviour can be anticipated when using the cell-periodic part of the eigenstates, see Eq. (8). In order to remark and clarify this issue, we have added a comment above Eq. (6) [recall that Eq. (8) is just an application of this] explaining how the wave-vector delta appears, as well as slightly enlarged the discussion below Eq. (12).

1.(c) Also, can any intuitive physical interpretation be given to Eqs. (11) and (12)? The latter looks like a current flow across the boundary, doesn't it? I suppose that for diagonal elements $\langle n\boldsymbol{k}|\hat{\boldsymbol{v}}|n\boldsymbol{k}\rangle_v$, the first term in (11) vanishes, right?

At a first sight, one surely thinks that the first term in Eq. (11) gives the contribution to velocity in the bulk while the second is the flow though the boundaries. However, recall that this is just a consequence of considering a finite volume and wavefunction. Therefore, we consider that such statement can be naive idea which may be far from a clear physical interpretation (contrarily as it happens in fluid mechanics, for instance), which is hard to do at the level of microscopic matrix elements. Gu et al. (where the equation is presented) do not mention anything about this detail, probably because no conclusions can be linked to this form of the equation. We have therefore avoided to make any comment that can potentially be a bit misleading. Regarding the second point, yes, as $\omega_{n\boldsymbol{k},n\boldsymbol{k}}=0$, all the contribution comes from the surface term. In this case one has that $\boldsymbol{C}_{nk,nk}=\partial_k \epsilon_n (\boldsymbol{k})$, accordingly to Eq. (9). We have included a discussion of this diagonal term in Sec. II.D.

1(d) I suspect Eqs. (11-12) can be recast in terms of integrals over a the interior and the boundary of a single primitive cell, with wave functions normalized to the unit cell, and with an average over $k$, where $k$ runs (as now) over the $N \times N \times N$ k-space mesh of the primitive Brillouin zone. It seems to me that this would be a more natural way to present the result.

One can try to reduce the integrals in the unit cell and see what happens. For instance, in a one dimensional model, the vertex of the integration line of length $L$ can be written $x_0=\eta a+\Delta$, where $\eta$ is an integer and $\Delta \in (-a/2,a/2)$. This serves to connect $x_0$ to a point inside the unit cell. We directly obtain \begin{equation} \begin{split} &\int_{x_0}^{x_0+L}dx\psi_{nk}^{(v)\ast}(x)x\psi_{n'k'}^{(v)}(x)= \ & \ \ \ \ \ e^{i(k'-k)\eta a}\int_{\Delta}^{\Delta+L}dx\psi_{nk}^{(v)\ast}(x)x\psi_{n'k'}^{(v)}(x) \ & \ \ \ \ \ +\eta a e^{i(k'-k)\eta a}\delta_{nn'}\delta_{kk'} \end{split} \end{equation} And similarly for the surface term. However, not much else can be done: as the integrand contains a non-periodic operator, we cannot reduce the integration over the whole space by a sum of integrals within a unit cell.

1.(e) Minor additional comment about Eq. (12): the final $\boldsymbol{r}$ is difficult to notice, but is crucial. Maybe there could be a way to call the reader'sattention to its presence.

We agree in this minor points and have added a sentence below Eq. (12) to remark it.

2. I found Eq. (7) and its discussion to be confusing. It seems to me that Eq. (7) is simply true without the last correction term. I think the authors have to find a way to be more explicit about what they really mean. For example, let $P_b$ be the projection operator onto the basis (``almost unity"), and define objects like $\hat{H}'=P_b \hat{H} P_b$, $\hat{p}'=P_b \hat{p} P_b$, $\hat{r}'=P_b \hat{r} P_b$, $\hat{v}'=P_b \hat{v} P_b$, etc. In this notation, I think what the authors are trying to say is that if the operators are replaced by their primed versions in Eq.(7), the $\Delta$ correction term is needed. Is that the idea? It needs to be clearer. By the way, I think there is a typo 7 lines below Eq. (7), where it should be Tr($r_\alpha, p_\beta)$=Tr($p_\beta,r_\alpha$) (when using primed operators in my notation).

We agree with the appreciation that such discussion needs to be clearer. As this effect has visible consequences in Figs. (2) and (3), we have kept the discussion but removed the delta term in the equation. The effect of the finite basis is discussed without any extra mathematical definitions for simplicity in the reading. We have also corrected the typo.

3. There is some confusing notation in Eq. (5). Putting $O(\boldsymbol{r})$ inside the integrand seems to suggest the $O$ is a local operator, but clearly this is not the case, or else $\hat{O}_{\boldsymbol{k}}$ would be the same as $\hat{O}$. I suppose what is meant would correspond to replacing

$$ \hat{O}(\boldsymbol{r})\psi_{n'k'} $$

by

$$ [\hat{O}(\boldsymbol{r})\psi_{n'k'}](\boldsymbol{r}), $$

or more pedantically, by

$$ \langle \boldsymbol{r}|\hat{O}|\psi_{n'k'}\rangle. $$

We understand the confusion and have rewritten Eq. (5) with the suggested notation.

4. Shortly after Eq. (5), the wording ¨such that still satisfies the Bloch theorem" is misleading. There was at least one more similar misuse later in the manuscript. Bloch's theorem is a theorem about the eigenstates of a Hamiltonian, which this wave function is not. The correct wording would instead be something like ¨of Bloch form". Even better would be to discuss whether $\hat{O}$ is a periodic operator, i.e., one that commutes with crystal translations; if so, then the product of such an operator with a wave function of Bloch form is automatically a wave function of Bloch form.

We also agree with this point and have used the term ``of Bloch form", which seems a natural way for us to describe those waves functions.

5. Another awkward detail of notation is the appearance of superscripts $(v)$ and $(n)$ in Eq. (19) and elsewhere. It looks as though the notation is parallel, whereas in fact the meanings of these superscripts are completely different.

We have made the change $c_{\alpha}^{(n)} (\boldsymbol{k}) \rightarrow c_{\alpha n}(\boldsymbol{k})$ to avoid posible confusions.

6. I think the LHS of the first Eq. (25) should be $\boldsymbol{A}_{\alpha \alpha'} (\boldsymbol{k})$. Also, the asymmetry between the two Eqs. (25) ¨looks wrong" although I think it is actually correct. It may be worth a few words of explanation. Perhaps it may be worth emphasizing that in the case of a nonorthogonal basis, A is not Hermitian, as it is in the orthonormal case. It also may be a good idea to explicitly write expression for the last term in Eq. (44) in a language parallel to the first term; I guess it comes to something like

$$ \sum_{\boldsymbol{R}}e^{i \boldsymbol{k} \cdot \boldsymbol{R}} \boldsymbol{R} \langle \alpha\boldsymbol{0}|\alpha' \boldsymbol{R}\rangle, $$

doesn't it?

We have added some commentary below Eq. (25) to highlight this points. Also, we agree and have done the second suggested change and identified this quantity with the gradient of the overlap matrix below Eq. (44).

7. The choice of the acronym CREN is not explained; what does it stand for?

To be honest: we do not really know. These basis set, originally introduced by Pacios and Christiansen (Ref. [15]), have stablished as one of the standard Gaussian-type basis sets used in Quantum Chemistry. They are usually labelled as CREN-BL or CREN-BS for large or small core pseudopotentials. For instance, the known database https://www.basissetexchange.org/ for GTO basis sets include the ``CREN-BL(BS)" basis for most elements. The information is found by clicking an element and looking in the list of available basis sets. We have included a reference to this online database page for completeness in our manuscript.