# A comprehensive study of the velocity, momentum and position matrix elements
for Bloch states using a local orbital basis

J. J. Esteve-Paredes  Departamento de Física de la Materia Condensada,
Universidad Autónoma de Madrid, E-28049 Madrid, Spain.  J. J. Palacios
Departamento de Física de la Materia Condensada, Universidad Autónoma de
Madrid, E-28049 Madrid, Spain.  Instituto Nicolás Cabrera (INC) and Condensed
Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid
Madrid, Spain.

###### Abstract

We present a comprehensive study of the velocity operator,
$\hat{\bm{v}}=\frac{i}{\hbar}[\hat{H},\hat{\bm{r}}]$, when used in crystalline
solids calculations. The velocity operator is key to the evaluation of a
number of physical properties and its computation, both from a practical and
fundamental perspective, has been a long-standing debate for decades. Our work
summarizes the different approaches found in the literature, connecting them
and filling the gaps in the sometimes non-rigorous derivations. In particular
we focus on the use of local orbital basis sets where the velocity operator
cannot be approximated by the $k$-derivative of the Bloch Hamiltonian matrix.
Among other things, we show how the correct expression can be found without
unequivocal mathematical steps, how the Berry connection makes its way in this
expression, and how to properly deal with the two popular gauge choices that
coexist in the literature. Finally, we explore its use in density functional
theory calculations by comparing with its real-space evaluation through the
identification with the canonical momentum operator. This comparison offers
us, in addition, a glimpse of the importance of non-local corrections, which
may invalidate the naive momentum-velocity correspondence.

††preprint: APS/123-QED

##  I Introduction

The quantum mechanical velocity $\hat{\bm{v}}$ operator plays a central role
in the evaluation of macroscopic optoelectronic properties of crystalline
solids. The velocity matrix elements (VMEs) are generically needed to
determine transitions between band states through several formalisms such as
Fermi’s golden rule Spataru  _et al._ (2005) for decay and optical excitation
processes or the more general Kubo linear response theory Mahan (1990).
Closely related to $\hat{\bm{v}}$, the canonical momentum operator
$\hat{\bm{p}}$ also plays a key role, but more from a methodological
standpoint. The momentum matrix elements (MME) are, for instance, needed to
find parameter-free effective models within $\bm{k}\cdot\hat{\bm{p}}$
perturbation theory Voon and Willatzen (2009).

It is common to consider the velocity operator as
$\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$ (in atomic units, which we will use
throughout the text), following a classical to quantum mechanics
identification through the Heisenberg equations of motion (Heisenberg
picture). In the cases where $\hat{H}$ only contains the kinetic energy and a
potential commuting with the position operator $\hat{\bm{r}}$ (i.e., in
absence of spin-orbit coupling or non-local potentials), to work with
$\hat{\bm{v}}$ or with $\hat{\bm{p}}$ becomes completely equivalent. In
coordinate representation this means that $i[H,\bm{r}]$ and
$-i\bm{\nabla}_{\bm{r}}$ are interchangable.

In the following, we review the state of the art of the uses and misuses of
these two operators as well as the position operator $\hat{\bm{r}}$ when
evaluating matrix elements between band states. Evaluating the commutator
matrix elements presents no problems when dealing with localized states, as in
atomic physics, but fundamental difficulties can be found when dealing with
Bloch eigenstates due to $\bm{r}\psi_{n\bm{k}}(\bm{r})$ not belonging to the
same Hilbert space as that of the states themselves, $\psi_{n\bm{k}}(\bm{r})$.
This issue has been addressed by Gu and coworkers Gu  _et al._ (2013), in
addition to presenting an extensive review of the existing $p$-$r$ relations
(as called in their work) in the literature. Gu et al. find the correct
relation between the momentum (or velocity for the case when $\hat{\bm{v}}$
and $\hat{\bm{p}}$ are equivalent) and position matrix elements:

$\begin{split}\bra{n\bm{k}}\hat{\bm{v}}\ket{n^{\prime}\bm{k}^{\prime}}_{v}=&\
i[\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k}^{\prime})]\bra{n\bm{k}}\hat{\bm{r}}\ket{n^{\prime}\bm{k}^{\prime}}_{v}\\\
&+\bm{C}_{n\bm{k},n^{\prime}\bm{k}^{\prime}}.\end{split}$ (1)

Here the matrix element is taken between eigenstates normalized to a finite
volume $v$ and obeying periodic boundary conditions (PBCs), therefore not
decaying at the boundaries even if the volume tends to infinity. The surface
term ${\bm{C}}$ is calculated on the surface of the solid. With this is mind,
one can easily convince oneself that the position (also called dipole) matrix
element depends on the origin of coordinates and that ${\bm{C}}$ compensates
this choice (as the VME cannot depend on the origin). Eq. (1) above presents a
fundamental view of $i[\hat{H},\hat{\bm{r}}]$ rather than convenient shortcut
to evaluate the VME, as a challenging integration in coordinate space is
needed on the right side of the equality. It also remarks the difference with
the atomic case, where the surface term does not appear.

In practice, Bloch eigenstates are very often represented in a Bloch basis
which, in turn, may be expanded in a local orbital basis. In this regard, a
good effort has been put in the actual evaluation of the VME over the last
decades. Between same-$k$ Bloch eigenstates, the VME can be calculated through
the following expression:

$\begin{split}\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=&\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\nabla_{\bm{k}}\braket{\alpha\bm{k}}{\hat{H}}{\alpha^{\prime}\bm{k}}_{v}\\\
&+\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\Big{[}i\epsilon_{n}(\bm{k})\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})\\\
&-i\epsilon_{n^{\prime}}(\bm{k})\bm{A}_{\alpha^{\prime}\alpha}^{\ast}(\bm{k})\Big{]},\end{split}$
(2)

where the $c$’s are the coefficients of the expansion of eigenstates in a
generic $\ket{\alpha\bm{k}}_{v}$ Bloch basis and
$\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})$ is the Berry connection associated
with such Bloch basis. The first term is sometimes referred to as the Peierls
approximation, while the second term is needed to deliver the full matrix
element. The suppression of the second term, as we will show later, can lead
to large and uncontrolled quantitative errors. This issue was first explored
by Pedersen et al. Pedersen  _et al._ (2001), by trying to complete the
Peierls approximation within a tight-binding scheme. We also point out the
remarkable work by Tomczak et al. Tomczak and Biermann (2009), introducing an
intra-unit cell correction to
$\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}(\bm{k})$ (we will see later that
this can be understood as a consequence of the gauge choice), and adding extra
terms to this quantity similar to those in Eq. (2). Tomczak _et al._
contribution was later replicated, perhaps in a clearer and more complete way,
in the work of Lee et al. Lee  _et al._ (2018), who presented a complete
expression similar to that in Eq. (2) for a general, nonorthonormal basis.
Actually, as originally reported in Ref. Lee  _et al._ (2018), the Berry
connection did not appear. The fact that Eq. (2) can be recasted in this form
will be shown below in this work, thus generalizing the evaluation of the VME
to any basis, not necessarily comprised of local orbitals.

In the light of Eq. (1), Eq. (2) presents a somewhat puzzling aspect: first,
there is no surface term and, second, no term depends on the placement of the
integration volume $v$. However, Eq. (2) was derived in Ref. Lee  _et al._
(2018) essentially in the same manner as Eq. (1) was derived in Ref. Gu  _et
al._ (2013), namely, by making use of $\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$
projected in a chosen basis, the coordinate basis in the former and a local
orbital basis in the latter. Addressing this intriguing observation is part of
our motivation to carry out the present work, as well as studying the role of
$\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}(\bm{k})$ in the calculation of
velocity matrix elements. We will also explain how the the popular expression
for position matrix elements given by Blount Blount (1962) fits into this
comparison.

To this end, we organize this study as follows. In Sec. II we present the main
theoretical ingredients by first recalling the differences between periodic
boundary conditions versus the infinite volume case when defining Bloch
eigenstates. We follow by introducing two ways of treating
$\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$, one relying on an integration in the
whole finite volume of the system, and the other based on using the $k$
representation for operators, involving matrix elements between the cell-
periodic part of the Bloch eigenstates. In Sec. III we explain how Eq. (2)
rigorously comes about from the second method, while showing the way it has
been previously derived in the literature is mathematically inconsistent, to
say the least. In Sec. IV we present a numerical study that gives us insight
into the quantitative error that one makes when assuming the equality
$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}=\braket{n\bm{k}}{\hat{\bm{p}}}{n^{\prime}\bm{k}^{\prime}}$,
even in the presence of a local potential, and the trade-off between
computational simplicity and accuracy when using the Peierls approximation in
a practical situation. Finally, we summarize our main conclusions in Sec. V.

##  II Theory

###  II.1 Preliminary definitions

We start by recalling the normalization choices for eigenstates in a crystal.
This turns out to be a key point to understand the relation between position
and velocity operators. Two options are compatible with Bloch theorem: one can
assume a finite volume $v$ normalization or let the eigenstates extend to all
space. In order to distinguish the two cases, we write

$\begin{split}\ket{n\bm{k}}_{v}=\frac{1}{\sqrt{N}}e^{i\bm{k}\cdot\hat{\bm{r}}}\ket{u_{n\bm{k}}},\\\
\ket{n\bm{k}}=\frac{1}{(2\pi)^{2/3}}e^{i\bm{k}\cdot\hat{\bm{r}}}\ket{u_{n\bm{k}}},\end{split}$
(3)

being the normalization conditions for every case case

$\begin{split}&\braket{n\bm{k}}{n^{\prime}\bm{k}^{\prime}}_{v}=\int_{v[\bm{x}_{0}]}d^{3}r\psi_{n\bm{k}}^{(v)\ast}(\bm{r})\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)\ast}(\bm{r})=\delta_{nn^{\prime}}\delta_{\bm{k}\bm{k}^{\prime}},\\\
&\braket{n\bm{k}}{n^{\prime}\bm{k}^{\prime}}=\int_{-\infty}^{\infty}d^{3}r\psi_{n\bm{k}}^{\ast}(\bm{r})\psi_{n^{\prime}\bm{k}^{\prime}}^{\ast}(\bm{r})=\delta_{nn^{\prime}}\delta(\bm{k}-\bm{k}^{\prime}).\end{split}$
(4)

Note that this convention assumes $\ket{u_{n\bm{k}}}$ to be normalized to one
in the finite volume case and to the unit cell volume (denoted with $\Omega$
in the following) in the distribution case. The $\int_{v[\bm{x}_{0}]}$ means
that the integration domain (the volume $v$) is defined by the parallelepiped
determined by $N_{i}\bm{a}_{i}$ with $N_{i}$ being the number of cells in each
direction given by the primitive vectors $\bm{a}_{i}$. Its origin vertex is
located at the $\bm{x}_{0}$ point, which have to be selected in the first
place. Finally $N=N_{1}N_{2}N_{3}$ is the total number of cells of the
crystal. Born-von Karman boundary conditions
$\psi_{n\bm{k}}(\bm{r}+N_{i}\bm{a}_{i})=\psi_{n\bm{k}}(\bm{r})$ are applied,
leading to a quantization of the crystal momentum according to
$\bm{k}=\frac{l_{1}}{N_{1}}\bm{G}_{1}+\frac{l_{2}}{N_{2}}\bm{G}_{2}+\frac{l_{3}}{N_{3}}\bm{G}_{3}$
with $\bm{G}$ being reciprocal lattice vectors and $l_{i}=-M_{i},\ldots,M_{i}$
such as $N_{i}=2M_{i}+1$. In general, matrix elements of operators whose
application on eigenstates breaks periodicity may depend on $\bm{x}_{0}$. On
the other hand, in the case of infinite volume normalization, the integrals
run over the whole unbounded space, including infinities. In this case,
$\bm{k}$ vectors become a dense set inside the Brillouin zone (BZ).

The representation of a given operator $\hat{O}$ in both cases becomes

$\begin{split}&O_{n\bm{k},n\bm{k}^{\prime}}^{(v)}\equiv\braket{n\bm{k}}{\hat{O}}{n^{\prime}\bm{k}^{\prime}}_{v}\\\
&\ \ \ \ \ \ \ \ \ \ \
=\int_{v[\bm{x}_{0}]}d^{3}r\psi_{n\bm{k}}^{(v)\ast}(\bm{r})O(\bm{r})\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r}),\\\
&O_{n\bm{k},n^{\prime}\bm{k}^{\prime}}\equiv\braket{n\bm{k}}{\hat{O}}{n^{\prime}\bm{k}^{\prime}}\\\
&\ \ \ \ \ \ \ \ \ \ \ =\int
d^{3}r\psi_{n\bm{k}}^{\ast}(\bm{r})O(\bm{r})\psi_{n^{\prime}\bm{k}^{\prime}}(\bm{r})\end{split}$
(5)

(we will assume that $\int$ is equivalent to $\int_{-\infty}^{\infty}$ in the
following). If the operator $\hat{O}$ is such that
$O(\bm{r})\psi_{n^{\prime}\bm{k}^{\prime}}(\bm{r})$ still satisfies Bloch
theorem, the matrix elements can be reduced to an integration within the unit
cell involving the periodic part of eigenstates

$\braket{n\bm{k}}{\hat{O}}{n^{\prime}\bm{k}^{\prime}}_{v}=\delta_{\bm{k}\bm{k}^{\prime}}\braket{u_{n\bm{k}}}{\hat{O}_{\bm{k}}}{u_{n^{\prime}\bm{k}}}_{\Omega}.$
(6)

where $\hat{O}_{\bm{k}}\equiv
e^{-i\bm{k}\cdot\hat{\bm{r}}}\hat{O}e^{i\bm{k}\cdot\hat{\bm{r}}}$, sometimes
called the “$k$ representation of an operator”. In the infinite volume case
the expression is similar but taking
$\delta_{\bm{k}\bm{k}^{\prime}}\rightarrow\delta(\bm{k}-\bm{k}^{\prime})/\Omega$.
In what follows we particularize to the velocity operator $\hat{\bm{v}}$ and
its relation to other quantities.

###  II.2 Relation between the velocity and momentum matrix elements

As discussed in the introduction, the velocity and momentum operators can only
be interchanged if spin-orbit coupling is neglected and the periodic potential
in the crystal is assumed to be local. Let the Hamiltonian be separated into
$\hat{H}=\hat{H}_{L}+\hat{H}^{\prime}$, where
$\hat{H}_{L}=\hat{\bm{p}}^{2}/2+V(\hat{\bm{r}})$ with $V(\hat{\bm{r}})$ the
local periodic part of the lattice potential and where no spin-orbit coupling
has been included. Then, for the first term one can write
$\hat{\bm{p}}=i[\hat{H}_{L},\hat{\bm{r}}]$, and the projection of the full
velocity operator on a subspace of band states can be written as

$\begin{split}\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}_{v}=&i\braket{n\bm{k}}{[\hat{H}_{L}+\hat{H}^{\prime},\hat{\bm{r}}]}{n^{\prime}\bm{k}^{\prime}}_{v}\\\
=&\braket{n\bm{k}}{\hat{\bm{p}}}{n^{\prime}\bm{k}^{\prime}}_{v}+i\braket{n\bm{k}}{[\hat{H}^{\prime},\hat{\bm{r}}]}{n^{\prime}\bm{k}^{\prime}}_{v}\\\
&+\Delta_{n\bm{k},n^{\prime}\bm{k}^{\prime}}^{(v)}.\end{split}$ (7)

$\Delta_{n\bm{k},n^{\prime}\bm{k}^{\prime}}^{(v)}$ is an extra term that one
encounters in practice and which is a consequence of the finite nature of the
Hilbert space in actual calculations. This can be ultimately traced back to
the fact that the canonical commutation relation
$[\hat{r}_{\alpha},\hat{p}_{\beta}]=i\delta_{\alpha\beta}$ can never exactly
hold for a finite-matrix representation since in such cases
$\text{Tr}(\hat{r}_{\alpha}\hat{p}_{\beta})=\text{Tr}(\hat{p}_{\alpha}\hat{r}_{\beta})$.
Therefore, one can only expect the quantity
$\Delta_{n\bm{k},n^{\prime}\bm{k}^{\prime}}^{(v)}$ to be negligible if the
physical states are sufficiently well represented in the working Hilbert
space. We will give below a few examples of this practical limitation.

The presence of the second term in Eq. (7) is challenging from a practical
standpoint and only when $\hat{H}^{\prime}=0$, Eq. (7) becomes the theoretical
velocity-momentum equality. In any case, evaluating the VME seems to require,
in principle, the evaluation of the MME through its representation
-$i\bm{\nabla}_{\bm{r}}$. In the following we explore two routes that can be
followed to by-pass the evaluation of the MME and, at the same time, of the
second term if present.

###  II.3 Relation between velocity matrix elements and the Berry connection

We first write the VME in the $k$ representation. It is easy to see that
$\hat{\bm{v}}_{\bm{k}}=e^{-i\bm{k}\cdot\hat{\bm{r}}}i[\hat{H},\hat{\bm{r}}]e^{i\bm{k}\cdot\hat{\bm{r}}}=\nabla_{\bm{k}}\hat{H}_{\bm{k}}$,
so one can write

$\bra{n\bm{k}}\hat{\bm{v}}\ket{n^{\prime}\bm{k}^{\prime}}_{v}=\delta_{\bm{k}\bm{k}^{\prime}}\braket{u_{n\bm{k}}}{(\nabla_{\bm{k}}\hat{H}_{\bm{k}})}{u_{n^{\prime}\bm{k}}}_{\Omega}.$
(8)

By applying chain rule it is straightforward to find

$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}_{v}=\delta_{\bm{k}\bm{k}^{\prime}}[i\omega_{n\bm{k},n^{\prime}\bm{k}}\bm{A}_{nn^{\prime}}(\bm{k})+\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}]$
(9)

where
$\omega_{n\bm{k},n^{\prime}\bm{k}}\equiv\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k})$
and $\bm{A}_{nn^{\prime}}(\bm{k})\equiv
i\braket{u_{n\bm{k}}}{\nabla_{\bm{k}}u_{n^{\prime}\bm{k}}}_{\Omega}$, this
last quantity being the Berry connection 111Note that Berry connection is not
an operator and hence it does not follow the initial definition of Eq. (5).
Eq. (9) can be found in the literature, see e.g. Ref. Ventura  _et al._
(2017). The equation above replaces Eq. (7) by introducing the evaluation of
the Berry connection associated with the Bloch eigenstates. This, however, can
be a cumbersome task since $k$ derivatives of eigenstates are not known in
numerical diagonalization procedures. While this problem can be circumvented
through perturbation theory Vanderbilt (2018), in Sec. III we show how Eq. (9)
can be recasted in a more convenient and familiar form. Incidentally, note
that Eq. (9) provides a way to compute the non-diagonal Berry connection
elements if the VMEs are known.

###  II.4 Relation between velocity and position matrix elements

Alternatively, we can directly perform the integrals that appear in Eq. (7)
when representing in coordinate space. Assuming $\hat{H}^{\prime}=0$ and,
therefore, being able to write $\hat{\bm{v}}\equiv
i[\hat{H},\hat{\bm{r}}]=\hat{\bm{p}}$, one is free to use
$-i\bm{\nabla}_{\hat{\bm{r}}}$ or $i[\hat{H},\hat{\bm{r}}]$. In both cases the
explicit knowledge of the real-space wavefunction of the eigenstates is
required. In the former case derivatives need to be carried out, which
depending on the orbital basis can be more or less cumbersome to implement. In
the latter, the use of the commutator entails further steps, where one needs
to pay attention to the correct use of the hermiticity of $\hat{H}$ in the
$\hat{H}\hat{\bm{r}}$ product. This procedure, which has been followed by Gu
and coworkers in Ref. Gu  _et al._ (2013), only applies to eigenstates in the
framework of finite volume normalization, where integrals for matrix elements
can be converged. One starts with

$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}_{v}=i\int_{v[\bm{x}_{0}]}d^{3}r\psi_{n\bm{k}}^{(v)\ast}(\bm{r})[H\bm{r}-\bm{r}H]\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r}),$
(10)

where we have to act with $-\nabla_{\bm{r}}^{2}$ on
$\bm{r}\psi_{n^{\prime}\bm{k}}^{(v)}(\bm{r})$ and
$\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r})$. After performing the
derivatives and using Gauss’s theorem, one obtains

$\begin{split}\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}_{v}=i\omega_{n\bm{k},n^{\prime}\bm{k}^{\prime}}\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}_{v}+\bm{C}_{n\bm{k},n^{\prime}\bm{k}^{\prime}},\end{split}$
(11)

where

$\begin{split}\bm{C}_{n\bm{k},n^{\prime}\bm{k}^{\prime}}=&-\frac{i}{2}\int_{\partial
v[\bm{x}_{0}]}d\bm{S}\cdot\text{\bigg{\\{}}\psi_{n\bm{k}}^{(v)\ast}(\bm{r})\nabla_{\bm{r}}\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r})\\\
&-[\nabla_{\bm{r}}\psi_{n\bm{k}}^{(v)}(\bm{r})]^{\ast}\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r})\text{\bigg{\\}}}\bm{r}\end{split}$
(12)

(same comment Note (1) applies here). The appearance of this last surface term
arises from the periodicity of the wave functions at the surface of the
material volume, which we denote with $\partial v[\bm{x}_{0}]$. It is
important to note that the wavefunctions do not decay even in the limit of an
infinite volume and this term is always present.

As noticed in Ref. Gu  _et al._ (2013), the hermiticity property cannot be
applied as usual in
$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}_{v}=i\braket{n\bm{k}}{[\hat{H}\hat{\bm{r}}-\hat{\bm{r}}\hat{H}]}{n^{\prime}\bm{k}^{\prime}}_{v}$,
which results in the surface term above. Secondly, both matrix elements on the
right hand side (RHS) in Eq. (11) depend on the origin of the integration
volume and are not $k$-diagonal, while the sum does not depend on this
arbitrary choice of origin and is diagonal in the wave vector as the VME
actually is. The relative weight of
$\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}_{v}$ versus
$\bm{C}_{n\bm{k},n^{\prime}\bm{k}^{\prime}}$ with respect to the full VME is
also explored in Ref. Gu  _et al._ (2013) showing that, in general, one cannot
find a point $\bm{x}_{0}$ that makes the surface term to vanish, even for
certain analytical models. Therefore, Eq. (11) presents no advantage versus
directly computing
$-i\braket{n\bm{k}}{\nabla_{\hat{\bm{r}}}}{n^{\prime}\bm{k}^{\prime}}$ to find
the VME, as one still has to perform nontrivial integrations for position and
surface matrix elements. It provides us, however, with the conclusion that
momentum and dipole matrix elements (multiplied by the frequency) should never
be interchanged when dealing with Bloch eigenstates in a finite volume.

###  II.5 Relation between velocity and position matrix elements with a
distribution basis

If Bloch eigenstates are normalized as distributions [recall Eqs. (3) and
(4)], then one can still use them as a basis to represent general physical
quantum states in the crystal. This was originally referred to as the crystal
momentum representation (CMR) Blount (1962), where one writes

$\ket{\phi}=\sum_{n}\int_{\text{BZ}}d^{3}kg_{n}(\bm{k})\ket{n\bm{k}},$ (13)

with $g_{n}(\bm{k})$ being the envelope function for the $n$ band. The matrix
elements between two physical states is written

$\braket{\phi_{1}}{\hat{O}}{\phi_{2}}=\sum_{nn^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}g_{n}^{(1)\ast}(\bm{k})g_{n^{\prime}}^{(2)}(\bm{k})\braket{n\bm{k}}{\hat{O}}{n^{\prime}\bm{k}^{\prime}}.$
(14)

Now one needs the find the matrix element
$\braket{n\bm{k}}{\hat{O}}{n^{\prime}\bm{k}^{\prime}}$ that enters the
calculation above. As only the full $n$ sums and $k$ integrations are
relevant, we can take into account the boundary properties of the state
$\ket{\phi}$. This is the case of the matrix elements for the position
operator, for which Blount Blount (1962) noticed that
$\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}$ is ill-defined by
itself but that a distribution form can be given if $\hat{\bm{r}}$ is assumed
to act on a state $\ket{\phi}$ belonging to its domain. Specifically, Blount
showed that

$\begin{split}\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}=&-i\nabla_{\bm{k}^{\prime}}\delta(\bm{k}^{\prime}-\bm{k})\delta_{nn^{\prime}}\\\
&+\frac{1}{\Omega}\delta(\bm{k}^{\prime}-\bm{k})\bm{A}_{nn^{\prime}}(\bm{k}).\end{split}$
(15)

The effect of boundary conditions on eigenstates is highlighted here, as
$\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}$ fundamentally
differs from $\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}_{v}$,
addressed in the previous section. Position operator matrix elements here do
not depend on any arbitrary origin, but only make sense within Eq. (14).

As far as the velocity operator is concerned, Eq. (9) is still perfectly valid
in the infinite volume case:

$\begin{split}\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}=&\delta(\bm{k}-\bm{k}^{\prime})[\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}\\\
&+\frac{i}{\Omega}\omega_{n\bm{k},n^{\prime}\bm{k}^{\prime}}\bm{A}_{nn^{\prime}}(\bm{k})],\end{split}$
(16)

expression to be used, again, only in the context of Eq. (14). Alternatively,
in Appendix A we also show that projecting $\hat{\bm{v}}$ on general physical
states
$\braket{\phi}{\hat{\bm{v}}}{\phi^{\prime}}=i\braket{\phi}{[\hat{H},\hat{\bm{r}}]}{\phi^{\prime}}$,
along with Eq. (15), also leads to Eq. (16).

Finally, to complete the connection between the different matrix element
expressions, it is straightforward to show that

$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}=i\omega_{n\bm{k},n^{\prime}\bm{k}^{\prime}}\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}^{\prime}}.$
(17)

Again, one simply needs to project
$\braket{\phi}{\hat{\bm{v}}}{\phi^{\prime}}=i\braket{\phi}{[\hat{H},\hat{\bm{r}}]}{\phi^{\prime}}$
and proceed in the same manner as explained in previous section. Here,
however, the surface term vanishes due to $\ket{\phi}$ and
$\ket{\phi^{\prime}}$ being square-integrable over all space. Notice that this
momentum and position relation matches that in atomic physics. This is also
shown in Ref. Gu  _et al._ (2013) by using narrow $k$ envelope functions in
the limit of zero width.

##  III Velocity matrix elements when representing in a Bloch basis.

Having established a comprehensive overview of the available recipes to
evaluate the VME and their proper use, we proceed now with their actual
computation when a generic and possibly non-orthonormal Bloch basis is used to
expand the Bloch eigenstates:

$\ket{n\bm{k}}_{v}=\sum_{\alpha}c_{\alpha\bm{k}}^{(n)}\ket{\alpha\bm{k}}_{v}.$
(18)

We stress again that $\ket{\alpha\bm{k}}_{v}$ is a generic basis state
satisfying Bloch’s theorem in a finite volume, with $\alpha$ being a generic
quantum number. The coefficients $c_{\alpha}^{(n)}(\bm{k})$ are found by
solving the generalized eigenvalue problem

$\sum_{\alpha^{\prime}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha^{\prime}}^{(n)}(\bm{k})=\epsilon_{n}(\bm{k})\sum_{\alpha^{\prime}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha^{\prime}}^{(n)}(\bm{k}),$
(19)

where $H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})$ and
$S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})$ are the matrices representing the
Hamiltonian and identity operators, respectively.

Eq. (9) can now be properly converted into more familiar expression. First the
Berry connection reads

$\begin{split}\bm{A}_{nn^{\prime}}(\bm{k})=&i\sum_{\alpha\alpha^{\prime}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha}^{(n)^{\ast}}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\\\
&+\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)^{\ast}}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}\bm{A}_{\alpha\alpha^{\prime}}(\bm{k}).\end{split}$
(20)

In this expression one has to perform derivatives in $k$ space of the
coefficients $c_{\alpha}^{(n)}(\bm{k})$. In most cases these coefficients are
obtained by numerical diagonalization of Eq. (19) so that they are not
continuous and, therefore, differentiable. However, this can be avoided by
directly employing the chain rule after inserting Eq. (20) into Eq. (9),
leading to

$\begin{split}&\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=\bm{v}_{nn^{\prime}}^{\text{(A,
cell)}}(\bm{k})+\bm{v}_{nn^{\prime}}^{\text{(B, cell)}}(\bm{k});\\\
&\bm{v}_{nn^{\prime}}^{\text{(A,
cell)}}(\bm{k})=\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}),\\\
&\bm{v}_{nn^{\prime}}^{\text{(B,
cell)}}(\bm{k})=i\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\Big{[}\epsilon_{n}(\bm{k})\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
-\epsilon_{n^{\prime}}(\bm{k})\bm{A}_{\alpha^{\prime}\alpha}^{\ast}(\bm{k})\Big{]}.\end{split}$
(21)

We show the complete derivation in Appendix B. Eq. (21) is one important
result of this work: it allows to compute VME from the Hamiltonian matrix
elements and the Berry connection in whichever Bloch basis. We have
differentiated two contributions, A and B, to the VME. The first one is evokes
the exact expression $\hat{\bm{v}}_{\bm{k}}=\nabla_{\bm{k}}\hat{H}_{\bm{k}}$,
but the second one is equally important, as we will show below. Eq. (21)
clearly shows that $\nabla_{\bm{k}}\hat{H}_{\bm{k}}$ is not, in general,
equivalent to $\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})$.

In many practical cases the Bloch basis is expanded, in turn, in a local
orbital basis. Regarding this, two different types of basis can be found in
the literature:

$\begin{split}\ket{\alpha\bm{k}}_{v}=&1/\sqrt{N}\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\ket{\alpha\bm{R}}_{v}\
\text{and}\\\
\ket{\tilde{\alpha}\bm{k}}_{v}=&1/\sqrt{N}\sum_{\bm{R}}e^{i\bm{k}\cdot(\bm{R}+\bm{d}_{\alpha})}\ket{\alpha\bm{R}}_{v},\end{split}$
(22)

where $\ket{\alpha\bm{R}}_{v}$ is an orbital with $\bm{d}_{\alpha}$ position
vector inside the unit cell of site $\bm{R}$. A finite size crystal containing
$N$ cells is assumed throughout. Bloch eigenstates are now given by

$\begin{split}&\ket{n\bm{k}}_{v}=\sum_{\alpha}c_{\alpha\bm{k}}^{(n)}\ket{\alpha\bm{k}}_{v},\
\text{or likewise}\\\
&\ket{n\bm{k}}_{v}=\sum_{\alpha}b_{\alpha\bm{k}}^{(n)}\ket{\tilde{\alpha}\bm{k}}_{v},\end{split}$
(23)

with both expansions being related by
$c_{\alpha\bm{k}}^{(n)}=e^{i\bm{k}\cdot\bm{d}_{\alpha}}b_{\alpha\bm{k}}^{(n)}$.
We will use the former basis in this work by default, which we refer to as the
cell gauge, and make considerations related to the other one, the atom gauge,
when appropriate. In this basis the matrices needed in Eq. (19) become

$\begin{split}H_{\alpha\alpha}^{(v)}(\bm{k})&=\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\hat{H}}{\alpha^{\prime}\bm{R}}\
\text{and}\\\
S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})&=\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\alpha^{\prime}\bm{R}}.\end{split}$
(24)

Now we can recast the Berry connection terms in Eq. (21) into an explicit form
involving the position operator

$\begin{split}&\bm{A}_{\alpha\alpha}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\hat{\bm{r}}}{\alpha^{\prime}\bm{R}}+i\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}),\\\
&\bm{A}_{\alpha^{\prime}\alpha}^{\ast}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\hat{\bm{r}}}{\alpha^{\prime}\bm{R}},\end{split}$
(25)

so that Eq. (21) becomes identical to that reported in Ref. Lee  _et al._
(2018) (see Appendix B for all the details.)

It is important to note that neither of the two terms in Eq. (21) is gauge
independent. One can easily check how the two terms change when switching to
the atom gauge, according to Eq. (22). For instance, the first term becomes

$\begin{split}\bm{v}_{nn^{\prime}}^{\text{(A,
atom)}}(\bm{k})&=\sum_{\alpha\alpha^{\prime}}b^{(n)\ast}_{\alpha}(\bm{k})b^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\nabla_{\bm{k}}\tilde{H}_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})\\\
&=\bm{v}_{nn^{\prime}}^{\text{(A, cell)}}(\bm{k})\\\ &\
+i\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})(\bm{d}_{\alpha^{\prime}}-\bm{d}_{\alpha}),\end{split}$
(26)

while the correction for the B term is the same with opposite sign, showing
that the absolute value of the sum
$\bm{v}_{nn^{\prime}}^{\text{(A)}}+\bm{v}_{nn^{\prime}}^{\text{(B)}}$ is gauge
invariant, as it should be for a physical operator.

It is worth obtaining the form of Eq. (21) in the maximally localized (Wannier
orbitals Marzari  _et al._ (2012)) or tight-binding limit, denoted here with
$\nu$. In this case, only the intra-atomic dipolar matrix elements in Eq. (25)
survive, leading to

$\begin{split}\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=&\sum_{\nu\nu^{\prime}}c_{\nu}^{(n)\ast}(\bm{k})c_{\nu^{\prime}}^{(n^{\prime})}(\bm{k})\nabla_{\bm{k}}H_{\nu\nu^{\prime}}(\bm{k})\\\
&+i\sum_{\nu\nu^{\prime}}c_{\nu}^{(n)\ast}(\bm{k})c_{\nu^{\prime}}^{(n^{\prime})}(\bm{k})H_{\nu\nu^{\prime}}(\bm{k})(\bm{d}_{\nu^{\prime}}-\bm{d}_{\nu}),\end{split}$
(27)

which is the same as Eq. (26). This tells us that, under a maximal
localization condition, computing the VME using only the gradient term (the A
term) in the atom gauge, therefore neglecting
$\bm{v}_{nn^{\prime}}^{(B,\text{atom})}$, is equivalent to computing both
terms (the full VME) in the cell gauge. This means that
$\bm{v}_{nn^{\prime}}^{(B,atom)}=0$, as can be easily checked. The second line
of Eq. (26) [or equivalently Eq. (27)] was presented in Ref. Wissgott  _et
al._ (2012) as a “Peierls substitution approach to the case of multiatomic
unit cells”. Based on our previous discussion, we see that it appears
naturally within the atom gauge. In the more general case of a non-orthonormal
basis, both terms of Eq. (21) must be evaluated regardless of the gauge
choice. We examine this more in depth in Sec. IV.

It is worth ending this section by briefly discussing the work of Lee et al.
Lee  _et al._ (2018). They present an expression for the VME which is, in
fact, a particular case of our general expression Eq. (21) (we reproduce it in
Appendix B). However, we believe that in order to reach their expression for
the VME, they have inadvertently mixed Hilbert spaces. Their derivation starts
from $\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$ and, briefly, they follow by
projecting
$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=i\braket{n\bm{k}}{[\hat{H},\hat{\bm{r}}]}{n^{\prime}\bm{k}}_{v}$,
expanding eigenstates in a non-orthonormal local orbital basis, and inserting
the closure relation
$\hat{I}=\sum_{\alpha\bm{R},\alpha^{\prime}\bm{R}^{\prime}}S_{\alpha\bm{R},\alpha^{\prime}\bm{R}^{\prime}}$
between the product of operators. We note that their procedure is equivalent
to start by writing

$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=i\omega_{n\bm{k},n^{\prime}\bm{k}}\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}},$
(28)

and proceeding in the same manner. The problem of starting with Eq. (28) is
that, as explained in Sec. II, this equality only holds in the case of open
boundary conditions (infinite systems). This is not the case when using local
orbital basis sets, where the band eigenstates obey a finite volume
normalization [see Eq. (5) and Eq. (22)]. The correct result found in Ref. Lee
_et al._ (2018) can only be explained by the unjustified identification of
$\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}}_{v}+\bm{C}_{n\bm{k},n^{\prime}\bm{k}^{\prime}}$
in Eq. (11) with
$\braket{n\bm{k}}{\hat{\bm{r}}}{n^{\prime}\bm{k}}=\int_{-\infty}^{\infty}d^{3}r\psi_{n\bm{k}}^{(v)\ast}(\bm{r})\bm{r}\psi_{n^{\prime}\bm{k}^{\prime}}^{(v)}(\bm{r})$,
which leads to Eq. (28). By doing this, the dependence on an arbitrary origin
of integration is effectively removed by the new integration limits, but the
integral is ill-defined. This subtle issue, which can be easily missed, is
stressed by our notation in Eq. (28), where we have put the subscript $v$ is
on the left hand side but not on the right hand side of the equality.

In the next section we present some numerical examples in order to explore the
details of the VME formulas in a practical situation.

##  IV Practical cases: hexagonal boron nitride and graphene

The first goal of this section is to gauge the importance of the different
terms in Eq. (7), by comparing between independent evaluations of the VME and
MME. Particularizing to a local orbital basis case, the former can be
evaluated from Eq. (21), while the latter becomes
$\braket{n\bm{k}}{\hat{\bm{p}}}{n^{\prime}\bm{k}}=-i\sum_{\alpha\alpha^{\prime}\bm{R}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\nabla_{\hat{\bm{r}}}}{\alpha^{\prime}\bm{R}}$.
Our second goal is to explore the relative importance of the two terms in Eq.
(21).

###  IV.1 Detailed numerical analysis of VME and MMe

We start by computing the band structure of a benchmark material. We choose a
monolayer of hexagonal boron nitride (hBN), which is a sufficiently complex
system to our purposes. In Fig. 1, we show: (i) a tight-binding (TB) two-band
calculation for the lower (upper) valence (conduction) bands, including only
first neighbour interactions between the $p_{z}$ orbitals of B and N atoms,
(ii) a DFT calculation employing a small-core pseudopotential basis set
Fernandez Pacios and Christiansen (1985) to replace the 1s2 electrons in every
atom (labelled as CREN) and, (iii) an all-electron calculation with the 6-31G*
basis set Dill and Pople (1975). The DFT calculations were performed using
CRYSTAL17 Dovesi  _et al._ (2018) with the local von Barth-Hedin exchange-
correlation functional von Barth and Hedin (1972). The use of Gaussian type
basis sets allows us to perform real-space integrations in an analytical
fashion. We are not concerned here with the accuracy of the obtained gap so we
have excluded the use of hybrid functionals and their possible extra non-local
contributions.

Both DFT band structures are essentially similar up to the conduction band.
The agreement is particularly good for the valence and conduction bands, both
with a band gap of 4.55 eV, except maybe for a noticeable difference at the M
point of $\simeq 0.5$ eV. As expected, only the more accurate all-electron
calculation with a large basis can reproduce results in the literature Galvani
_et al._ (2016). The tight-binding parameters can be fitted to resemble one of
these calculations. It is easier to obtain a better overall fit to the CREN
band structure with $t=2.15$ eV, as only two $p_{z}$ orbitals are present to
reproduce the energy dispersion.

Figure 1: Comparison of the band structure of monolayer hBN using different
approaches: (i) a first-neighbour tight-binding two-band model with 2.15 eV
hopping, (ii) a DFT calculation using a small-core pseudopotential basis set
and (iii) a DFT all-electron calculation (see text for further details).

We now explore in some detail Eqs. (7) and (21). To this purpose, in Fig. 2 we
show the magnitude of several quantities relevant to the band-gap optical
transition along the $\Gamma-\text{K}-\text{M}$ path. Fig. 2(a) shows the
absolute value of the $x$ component of the VME and the MME for the three cases
shown in Fig. 1. Looking at the VME, the TB result deviates quantitatively
from the other two, but not qualitatively. When comparing the VME and the MME,
we observe that for the CREN basis the difference is significant, particularly
near M, while that for the large basis this difference is negligible. The
latter result proves that Eq. (21) is properly implemented since the
evaluation of the MME is essentially analytical due to the use of Gaussian
orbitals. The difference found in the former calculation can be attributed, as
reflected in Eq. (7), to both the presence of the non-local pseudopotential
and the difference in the size of the Hilbert space (8 bands versus 36 bands).
In order to isolate the effect of the non-locality from that of the basis
size, we have repeated the calculation with the all-electron basis using the
non-local functional HSE06 Krukau  _et al._ (2006) (not shown in the figures).
In this case, the deviation between the VME and MME curves becomes
appreciable, but not larger than the one for the CREN case. In summary, these
results explicitly show that the VME and MME cannot always be taken as the
same quantity.

In Fig. 2(b), we compare the magnitude of the $k$-gradient term [the A term in
Eq. (21)] calculated in the atom gauge for the three different cases. In the
TB case, this term gives the full value for the VME. In DFT the results
deviate significantly from the exact value, showing the importance of the B
term in Eq. (21). The CREN basis presents a larger deviation, despite being
smaller in size [see Eq. (25)]. As mentioned in Sec. III, only the maximal
localization condition for the basis orbitals ensures that $v^{(A,\text{
atom})}$ gives the exact VME. This condition is not met in neither of the two
DFT basis sets used in our calculations. We also show the result obtained in
the cell gauge in Fig. 2(c). Now, not only quantitative differences appear,
but also selection rules break when approaching the $\Gamma$ point (here the
VME must be zero according to the irreducible representations of the wave
functions). Therefore, identifying the VME simply as a $k$-derivative of the
Bloch Hamiltonian in the cell gauge can lead, not only to quantitative errors,
but also to incorrect physical interpretations.

Figure 2: Absolute value of matrix elements for the band-gap transition along
the Brillouin zone of monolayer hBN. (a) Velocity and momentum matrix elements
for the two DFT calculations presented in 1 and in the tight-binding
approximation (b) Same for the first term of Eq. (21) in the atom gauge [see.
Eq. (23)]. (c) Same as (b) but in the cell gauge.

Figure 3: Frequency-dependent sheet conductivity of monolayer hBN as obtained
from the evaluation of Eq. (29) using (a) a all-electron and (b) small-core
pseudopotential DFT calculations. The calculation with a first-neighbour
tight-binding model is included in both panels. See Fig. 1 for the
corresponding band structures.

###  IV.2 Optical conductivity

The calculation of an experimentally measurable quantity such as the optical
conductivity can be affected by an incorrect evaluation of the VME. To show
this we make use of the Kubo-Greenwood Louie and Cohen (2006) expression (we
do not use atomic units here for clarity):

$\begin{split}\sigma_{\alpha\beta}(\omega)=&-\frac{i}{N_{k}}\frac{e^{2}\hbar}{\Omega}\sum_{nn^{\prime}\bm{k}}\left(\frac{f_{n\bm{k}}-f_{n^{\prime}\bm{k}}}{\omega_{n\bm{k},n^{\prime}\bm{k}}}\right)\\\
&\times\frac{\braket{n\bm{k}}{\hat{v}_{\alpha}}{n^{\prime}\bm{k}}_{v}\braket{n^{\prime}\bm{k}}{\hat{v}_{\beta}}{n\bm{k}}_{v}}{\hbar\omega+\omega_{n\bm{k},n^{\prime}\bm{k}}+i\hbar\eta},\end{split}$
(29)

where $f_{n\bm{k}}$ is the Fermi-distribution occupation number and $N_{k}$ is
the number of $k$ points in the discretized Brillouin zone.

In Fig. 3 we show the longitudinal optical conductivity, computing the VMEs
within the different approximations considered in previous section. We have
separated the results obtained with the small-core basis from those with the
all-electron basis, as shown in Fig. 3(a) and Fig. 3(b), respectively, where
we have also added the calculation with the TB model in both panels. At the
bandgap frequency, the DFT and TB calculations involving the exact VME are
able to reproduce the quasiuniversal behaviour Stauber  _et al._ (2015) for a
parabolic noninteracting semiconductor, yielding $\sigma=e^{2}/2\hbar$. The
use of MMEs, instead of the VMES, fails for the pseudopotential and small
basis case [black dashed line in Fig. 3(b)], as expected from the discussion
in previous subsection. At higher frequencies the TB model underestimates the
response, which is similar in magnitude for both DFT cases, the only
difference being the position of the Van Hove singularity which originates in
the bands at the M point (see Fig. 1). In both DFT calculations the
$k$-gradient approximation overestimates the exact result. A calculation with
the $k$-gradient term in the cell gauge
$\hat{\bm{v}}\rightarrow\bm{v}^{(A,\text{ cell})}$ (not shown) gives an even
larger discrepancy at all frequencies, but worse, also removes the isotropic
behaviour of the conductivity tensor with $\sigma_{xx}\neq\sigma_{yy}$. This
erroneous behaviour has been already discussed for graphene in Ref. Nguyen and
Nguyen (2016) and highlights the importance of taking the $k$-derivative
approximation for VME using the appropriate gauge. It is also worth mentioning
here the work by Wissgott et al. Wissgott  _et al._ (2012). There, the Peierls
approximation in the atom gauge is tested versus the complete VME also through
a conductivity analysis of transition-metal oxides. Our conclusion about the
gauge choice, not explored in their work, could give a better insight about
the discrepancies that are found in Ref. Wissgott  _et al._ (2012).

A direct comparison with experiments can be made by analyzing the optical
response of graphene. It is known that monolayer graphene shows a quasi-
constant absorbance of $\sim 2.3\%$, corresponding to $\sigma=e^{2}/4\hbar$,
over the energy region that goes from the far-infrared to the visible spectrum
where excitonic effects are negligible ($<2$ eV) Mak  _et al._ (2008, 2011);
Yang  _et al._ (2009). Therefore, in this energy range, Kubo-Greenwood DFT-
based calculations are expected to give a faithful optical response. In Fig. 4
we show the optical conductivity calculated with two different basis sets,
equivalent to those used for hBN. We present results for the exact VMEs and
their approximated values using the MMEs. Experimental results from Ref. Mak
_et al._ (2008) are also shown. We can see that both basis sets give results
in very good agreement with the experimental ones when employing VME. For the
case of MMEs, the CREN basis set gives $\sim 0.175e^{2}/\hbar$, which
translates in a 30 % error when comparing to the experimental curve. This
result complements our previous study of hBN, showing the significant effect
of non-local operators and finite basis sets when trying to replace the VMEs
by the MMEs.

Figure 4: Same as Fig. 3 for the case of graphene. Experimentals results from
Ref. Mak  _et al._ (2008) are shown. VME and MMe has been used to represent
the velocity operator in the Kubo-Greenwood formula for the two DFT
calculations (other approximations are not shown in this case).

##  V Conclusion

We have presented a comprehensive study of the evaluation of VME in
crystalline solids, as obtained from the fundamental relation
$\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$. We have scrutinized all available
expressions in the literature, filling the gaps and connecting them in a
coherent story. We have seen that, when working in coordinate representation,
one is bound to deal with a very inconvenient surface term which can be
avoided by going first into the $k$-representation. We have obtained a general
expression which contains a familiar $k$-derivative term plus a correction
term which involves the Berry connection of the Bloch basis elements. When
using local orbitals as a basis, this can be rewritten in a more familiar form
(see, e.g., Ref. Lee  _et al._ (2018)), but whose previous derivations contain
unjustified mathematical steps. We have also shown several equivalences which
involve the momentum and position operators, including well-known expressions
in the crystal momentum representation (nonphysical distribution basis).

We have numerically tested the validity of different approximations to the VME
by computing the optical conductivity of monolayer hBN and graphene through
the Kubo-Greenwood formula. In particular, we have shown that approximating
the VME by $\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}(\bm{k})$ in a non-
orthonormal basis produces significant quantitative errors and may also give
rise to qualitative ones if one is not careful with the choice of gauge. We
have also made emphasis on the fact that the velocity and momentum matrix
elements can only be safely interchanged if the Hamiltonian is free of non-
local terms and eigenstates are well-represented in the working Hilbert space.

In summary, this work may well serve as a complete as well as rigorous guide
to the intricate relations behind the evaluation of the velocity, momentum,
and position matrix elements in crystalline solids.

###### Acknowledgements.

The authors acknowledge financial support from Spanish MINECO through Grant
No. PID2019- 109539GB-C43, the María de Maeztu Program for Units of Excellence
in R&D (Grant No. CEX2018-000805-M), the Comunidad Autónoma de Madrid through
the Nanomag COST-CM Program (Grant No. S2018/NMT-4321), the Generalitat
Valenciana through Programa Prometeo/2021/01, the Centro de Computación
Científica of the Universidad Autónoma de Madrid and the computer resources of
the Red Española de Supercomputación.

##  Appendix A Representation of $\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$ in the
distribution basis

In this appendix we prove that projecting
$\hat{\bm{v}}=i[\hat{H},\hat{\bm{r}}]$ in CMR along with the corresponding
expression for the position operator, Eq. (15), allows to obtain Eq. (9). Let
$\ket{m\bm{k}}$ be a general Bloch basis (orthonormal for simplicity)
following a distribution normalization. An identification with the eigenstates
basis will be made in the end. We have

$\braket{\phi_{1}}{\hat{\bm{v}}}{\phi_{2}}=i\braket{\phi_{1}}{[\hat{H},\hat{\bm{r}}]}{\phi_{2}}=i\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})\braket{m\bm{k}}{[\hat{H}\hat{\bm{r}}-\hat{\bm{r}}\hat{H}]}{m^{\prime}\bm{k}^{\prime}}.$
(30)

Now we insert the closure relation between the two operators:

$\begin{split}\braket{\phi_{1}}{\hat{\bm{v}}}{\phi_{2}}=&i\braket{\phi_{1}}{[\hat{H},\hat{\bm{r}}]}{\phi_{2}}\\\
=&i\sum_{mm^{\prime}m^{\prime\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}d^{3}k^{\prime\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})\times(H_{m\bm{k},m^{\prime\prime}\bm{k}^{\prime\prime}}\bm{r}_{m^{\prime\prime}\bm{k}^{\prime\prime},m^{\prime}\bm{k}^{\prime}}-\bm{r}_{m\bm{k},m^{\prime\prime}\bm{k}^{\prime\prime}}H_{m^{\prime\prime}\bm{k}^{\prime\prime},m^{\prime}\bm{k}^{\prime}})\\\
\equiv&\braket{\phi_{1}}{\hat{\bm{v}}^{(1)}}{\phi_{2}}+\braket{\phi_{1}}{\hat{\bm{v}}^{(2)}}{\phi_{2}}.\end{split}$
(31)

We have splitted the full matrix elements into two terms according to the two
parts in Eq. (15). First we work out
$\braket{\phi_{1}}{\hat{\bm{v}}^{(1)}}{\phi_{2}}$,

$\begin{split}\braket{\phi_{1}}{\hat{\bm{v}}^{(1)}}{\phi_{2}}=&\frac{i}{\Omega}\sum_{mm^{\prime}m^{\prime\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}d^{3}k^{\prime\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})H_{mm^{\prime\prime}}(\bm{k})\delta(\bm{k}-\bm{k}^{\prime\prime})[-i\delta_{m^{\prime\prime}m^{\prime}}\nabla_{\bm{k}^{\prime}}\delta(\bm{k}^{\prime}-\bm{k}^{\prime\prime})]\\\
&-\frac{i}{\Omega}\sum_{mm^{\prime}m^{\prime\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}d^{3}k^{\prime\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})[-i\delta_{mm^{\prime\prime}}\nabla_{\bm{k}^{\prime\prime}}\delta(\bm{k}^{\prime\prime}-\bm{k})]H_{m^{\prime\prime}m^{\prime}}(\bm{k}^{\prime\prime})\delta(\bm{k}^{\prime\prime}-\bm{k}^{\prime}).\end{split}$
(32)

We have taken into account that a crystal Hamiltonian is diagonal in the
$\bm{k}$ vector, this is
$H_{m\bm{k},m^{\prime}\bm{k}^{\prime}}\equiv\Omega^{-1}H_{mm^{\prime}}(\bm{k})\delta(\bm{k}-\bm{k}^{\prime})$.
Using the identity
$F(\bm{k})[\nabla_{\bm{k}}\delta(\bm{k}-\bm{k}^{\prime})]=-[\nabla_{\bm{k}}F(\bm{k})]\delta(\bm{k}-\bm{k}^{\prime})$
straightforwardly, one can see

$\begin{split}\braket{\phi_{1}}{\hat{\bm{v}}^{(1)}}{\phi_{2}}=&\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}d^{3}k^{\prime\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})H_{mm^{\prime}}(\bm{k})\delta(\bm{k}-\bm{k}^{\prime\prime})[\nabla_{\bm{k}^{\prime}}\delta(\bm{k}^{\prime}-\bm{k}^{\prime\prime})]\\\
&-\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}d^{3}k^{\prime\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})[\nabla_{\bm{k}^{\prime\prime}}\delta(\bm{k}^{\prime\prime}-\bm{k})]H_{mm^{\prime}}(\bm{k}^{\prime\prime})\delta(\bm{k}^{\prime\prime}-\bm{k}^{\prime})\\\
=&-\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kg_{m}^{(1)\ast}(\bm{k})[\nabla_{\bm{k}}g_{m^{\prime}}^{(2)}(\bm{k})]H_{mm^{\prime}}(\bm{k})+\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kg_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k})[\nabla_{\bm{k}}H_{mm^{\prime}}(\bm{k})]\\\
&-\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}\nabla_{\bm{k}}[g_{m}^{(1)\ast}(\bm{k})H_{mm^{\prime}}(\bm{k})]g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})\delta(\bm{k}-\bm{k}^{\prime}),\end{split}$
(33)

and applying the chain rule,

$\braket{\phi_{1}}{\hat{\bm{v}}^{(1)}}{\phi_{2}}=-\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}k\nabla_{\bm{k}}[g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k})H_{mm^{\prime}}(\bm{k})]+\frac{1}{\Omega}\sum_{mm^{\prime}}\int_{\text{BZ}}d^{3}kg_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k})[\nabla_{\bm{k}}H_{mm^{\prime}}(\bm{k})].$
(34)

The first term is zero following the conditions required by Blount Blount
(1962). Now we look at the Berry connection term

$\begin{split}\braket{\phi_{1}}{\hat{\bm{v}}^{(2)}}{\phi_{2}}&=\frac{i}{\Omega^{2}}\sum_{mm^{\prime}m^{\prime\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}g_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})\\\
&\times[H_{mm^{\prime\prime}}(\bm{k})\bm{A}_{m^{\prime\prime}m^{\prime}}(\bm{k}^{\prime\prime})\delta(\bm{k}-\bm{k}^{\prime\prime})\delta(\bm{k}^{\prime}-\bm{k}^{\prime\prime})-\bm{A}_{mm^{\prime\prime}}(\bm{k})H_{m^{\prime\prime}m^{\prime}}(\bm{k}^{\prime\prime})\delta(\bm{k}^{\prime\prime}-\bm{k})\delta(\bm{k}^{\prime\prime}-\bm{k}^{\prime})]\\\
&=\frac{i}{\Omega^{2}}\sum_{mm^{\prime}m^{\prime\prime}}\int_{\text{BZ}}d^{3}kg_{m}^{(1)\ast}(\bm{k})g_{m^{\prime}}^{(2)}(\bm{k}^{\prime})[H_{mm^{\prime\prime}}(\bm{k})\bm{A}_{m^{\prime\prime}m^{\prime}}(\bm{k})-\bm{A}_{mm^{\prime\prime}}(\bm{k})H_{m^{\prime\prime}m^{\prime}}(\bm{k})].\end{split}$
(35)

We now find the expression in the eigenstates basis. For clarity we rename
$m=n$, and use
$H_{nn^{\prime}}(\bm{k})=\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}\Omega$,
obtaining

$\begin{split}\braket{\phi_{1}}{\hat{\bm{v}}}{\phi_{2}}&=\sum_{nn^{\prime}}\int_{\text{BZ}}d^{3}kg_{n}^{(1)\ast}(\bm{k})g_{n^{\prime}}^{(2)}(\bm{k})\Big{\\{}\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}+\frac{i}{\Omega}[\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k})]\bm{A}_{nn^{\prime}}(\bm{k})\Big{\\}}\\\
&\equiv\sum_{nn^{\prime}}\int_{\text{BZ}}d^{3}kd^{3}k^{\prime}g_{n}^{(1)\ast}(\bm{k})g_{n^{\prime}}^{(2)}(\bm{k}^{\prime})\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}},\end{split}$
(36)

where
$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}^{\prime}}=\delta(\bm{k}-\bm{k}^{\prime})\left[\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}+i\Omega^{-1}\omega_{n\bm{k},n^{\prime}\bm{k}^{\prime}}\bm{A}_{nn^{\prime}}(\bm{k})\right]$,
which is precisely Eq. (16).

##  Appendix B Derivation of Eq. (21)

We start from Eq. (9) for the case $\bm{k}=\bm{k}^{\prime}$,

$\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=i\omega_{n\bm{k},n^{\prime}\bm{k}}\bm{A}_{nn^{\prime}}(\bm{k})+\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}$
(37)

Recall that the Berry connection is defined
$\bm{A}_{nn^{\prime}}(\bm{k})\equiv
i\braket{u_{n\bm{k}}}{\nabla_{\bm{k}}u_{n^{\prime}\bm{k}}}_{\Omega}$.
Expanding the periodic part in a Bloch Basis, $\ket{n\bm{k}}$ state is
$\ket{u_{n\bm{k}}}=\sum_{\alpha}c_{\alpha}^{(n)}(\bm{k})\ket{u_{\alpha\bm{k}}}$,
we readily obtain

$\begin{split}\bm{A}_{nn^{\prime}}(\bm{k})=&i\sum_{\alpha\alpha^{\prime}}S_{\alpha\alpha^{\prime}}(\bm{k})c_{\alpha}^{(n)^{\ast}}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\\\
&+\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)^{\ast}}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}\bm{A}_{\alpha\alpha^{\prime}}(\bm{k}).\end{split}$
(38)

We now insert this expression into Eq. (37), obtaining

$\begin{split}&\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=\\\ &\ \
-\epsilon_{n}(\bm{k})\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})\\\
&\ \
+\epsilon_{n^{\prime}}(\bm{k})\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})+\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}\\\
&\ \
+i[\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k})]\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\bm{A}_{\alpha\alpha^{\prime}}(\bm{k}).\end{split}$
(39)

Applying the chain rule in the second term

$\begin{split}&\epsilon_{n^{\prime}}(\bm{k})\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})=\\\
&\ \ \
-\epsilon_{n^{\prime}}(\bm{k})\sum_{\alpha\alpha^{\prime}}\nabla_{\bm{k}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})\\\
&\ \ \
-\epsilon_{n^{\prime}}(\bm{k})\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha}^{(n^{\prime})}(\bm{k})\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}),\end{split}$
(40)

so we have

$\begin{split}&\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=\\\ &\ \
-\sum_{\alpha\alpha^{\prime}}\epsilon_{n}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha}^{(n)\ast}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\\\
&\ \
-\sum_{\alpha\alpha^{\prime}}\nabla_{\bm{k}}c_{\alpha}^{(n)}(\bm{k})\epsilon_{n^{\prime}}(\bm{k})S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})\ast}(\bm{k})\\\
&\ \
-\epsilon_{n^{\prime}}(\bm{k})\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})+\nabla_{\bm{k}}\epsilon_{n}(\bm{k})\delta_{nn^{\prime}}\\\
&\ \
+i[\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k})]\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\bm{A}_{\alpha\alpha^{\prime}}(\bm{k}).\end{split}$
(41)

Now we can introduce the Hamiltonian matrix elements in the first two terms
according to the eigenvalue equation, yielding

$\begin{split}&-\sum_{\alpha\alpha^{\prime}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha}^{(n)\ast}(\bm{k})\nabla_{\bm{k}}c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\\\
&-\sum_{\alpha\alpha^{\prime}}\nabla_{\bm{k}}c_{\alpha}^{(n)\ast}(\bm{k})H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})=\\\
&\ \ \
-\nabla_{\bm{k}}\Bigg{[}\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\Bigg{]}\\\
&\ \ \
+\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}).\end{split}$
(42)

The first term cancels the gradient of the energy band in Eq. (41). In order
to write the final form of the expression, we note that
$\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})\equiv
i\braket{u_{\alpha\bm{k}}}{\nabla_{\bm{k}}u_{\alpha^{\prime}\bm{k}}}_{\Omega}=i\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k})+\bm{A}_{\alpha^{\prime}\alpha}^{\ast}(\bm{k})$,
which leave us with

$\begin{split}&\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=\bm{v}_{nn^{\prime}}^{\text{(A)}}(\bm{k})+\bm{v}_{nn^{\prime}}^{\text{(B)}}(\bm{k});\\\
&\bm{v}_{nn^{\prime}}^{\text{(A)}}(\bm{k})=\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}),\\\
&\bm{v}_{nn^{\prime}}^{\text{(B)}}(\bm{k})=\sum_{\alpha\alpha^{\prime}}c^{(n)\ast}_{\alpha}(\bm{k})c^{(n^{\prime})}_{\alpha^{\prime}}(\bm{k})\Big{[}i\epsilon_{n}(\bm{k})\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})\\\
&\ \ \ \ \ \ \ \ \ \ \
-i\epsilon_{n^{\prime}}(\bm{k})\bm{A}_{\alpha^{\prime}\alpha}^{\ast}(\bm{k})\Big{]},\end{split}$
(43)

as presented in the main text. Finally this expression is particularized for
Bloch states expanded in a local orbital basis, where
$\ket{\alpha\bm{k}}_{v}=1/\sqrt{N}\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\ket{\alpha\bm{R}}_{v}$,
leading to the Berry connection

$\bm{A}_{\alpha\alpha^{\prime}}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\hat{\bm{r}}}{\alpha^{\prime}\bm{R}}+i\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}^{(v)}(\bm{k}).$
(44)

The expression above is a generalization for that of a Wannier basis, see e.g.
Ref. Vanderbilt (2018); Marzari  _et al._ (2012). Here, an extra term arises
accounting from the nonorthonormal character of atomic states, differently
from the Wannier orbitals, which are orthonormal by construction. This is also
reflected by the appearance of the overlap matrix in the first line of Eq.
(20). In Eq. (44), dipole matrix elements between the basis set are integrated
in all space and not in the unit cell, different than in the original
definition for the Berry connection. This change is done by passing from
$\int_{\text{cell}}$ to $\lim_{N\to\infty}\frac{1}{N}\int_{-\infty}^{\infty}$,
that is well-defined for a periodic integrand. Finally Eq. (43) can be written

$\begin{split}&\braket{n\bm{k}}{\hat{\bm{v}}}{n^{\prime}\bm{k}}_{v}=\\\ &\ \ \
\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha^{\prime}}^{(n^{\prime})}(\bm{k})\Big{[}\nabla_{\bm{k}}H_{\alpha\alpha^{\prime}}(\bm{k})-\epsilon_{n}(\bm{k})\nabla_{\bm{k}}S_{\alpha\alpha^{\prime}}(\bm{k})\Big{]}\\\
&\ \ \
+i[\epsilon_{n}(\bm{k})-\epsilon_{n^{\prime}}(\bm{k})]\sum_{\alpha\alpha^{\prime}}c_{\alpha}^{(n)\ast}(\bm{k})c_{\alpha}^{(n^{\prime})}(\bm{k})\\\
&\ \ \
\times\sum_{\bm{R}}e^{i\bm{k}\cdot\bm{R}}\braket{\alpha\bm{0}}{\hat{\bm{r}}}{\alpha^{\prime}\bm{R}},\end{split}$
(45)

which is the formula given in Ref. Lee  _et al._ (2018).

## References

  * Spataru  _et al._ (2005) C. D. Spataru, S. Ismail-Beigi, R. B. Capaz, and S. G. Louie, Phys. Rev. Lett. 95, 247402 (2005).
  * Mahan (1990) G. Mahan,  _Many-Particle Physics_, Physics of Solids and Liquids (Springer US, 1990).
  * Voon and Willatzen (2009) L. Voon and M. Willatzen,  _The k p Method: Electronic Properties of Semiconductors_ (Springer Berlin Heidelberg, 2009).
  * Gu  _et al._ (2013) B. Gu, N. H. Kwong, and R. Binder, Phys. Rev. B 87, 125301 (2013).
  * Pedersen  _et al._ (2001) T. G. Pedersen, K. Pedersen, and T. Brun Kriestensen, Phys. Rev. B 63, 201101 (2001).
  * Tomczak and Biermann (2009) J. M. Tomczak and S. Biermann, Phys. Rev. B 80, 085117 (2009).
  * Lee  _et al._ (2018) C.-C. Lee, Y.-T. Lee, M. Fukuda, and T. Ozaki, Phys. Rev. B 98, 115115 (2018).
  * Blount (1962) E. I. Blount, Solid State Phys. 13, 305 (1962).
  * Note (1) Note that Berry connection is not an operator and hence it does not follow the initial definition of Eq. (5).
  * Ventura  _et al._ (2017) G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes, and N. M. R. Peres, Phys. Rev. B 96, 035431 (2017).
  * Vanderbilt (2018) D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Wissgott  _et al._ (2012) P. Wissgott, J. Kuneš, A. Toschi, and K. Held, Phys. Rev. B 85, 205133 (2012).
  * Fernandez Pacios and Christiansen (1985) L. Fernandez Pacios and P. A. Christiansen, The Journal of Chemical Physics 82, 2664 (1985).
  * Dill and Pople (1975) J. D. Dill and J. A. Pople, The Journal of Chemical Physics 62, 2921 (1975).
  * Dovesi  _et al._ (2018) R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro, and B. Kirtman, WIREs Computational Molecular Science 8, e1360 (2018).
  * von Barth and Hedin (1972) U. von Barth and L. Hedin, Journal of Physics C: Solid State Physics 5, 1629 (1972).
  * Galvani  _et al._ (2016) T. Galvani, F. Paleari, H. P. C. Miranda, A. Molina-Sánchez, L. Wirtz, S. Latil, H. Amara, and F. m. c. Ducastelle, Phys. Rev. B 94, 125303 (2016).
  * Krukau  _et al._ (2006) A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, The Journal of Chemical Physics 125, 224106 (2006).
  * Louie and Cohen (2006) S. Louie and M. Cohen,  _Conceptual Foundations of Materials: A Standard Model for Ground- and Excited-State Properties_, ISSN (Elsevier Science, 2006).
  * Stauber  _et al._ (2015) T. Stauber, D. Noriega-Pérez, and J. Schliemann, Phys. Rev. B 91, 115407 (2015).
  * Nguyen and Nguyen (2016) H.-V. Nguyen and V. H. Nguyen, Phys. Rev. B 94, 117401 (2016).
  * Mak  _et al._ (2008) K. F. Mak, M. Y. Sfeir, Y. Wu, C. H. Lui, J. A. Misewich, and T. F. Heinz, Phys. Rev. Lett. 101, 196405 (2008).
  * Mak  _et al._ (2011) K. F. Mak, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 106, 046401 (2011).
  * Yang  _et al._ (2009) L. Yang, J. Deslippe, C.-H. Park, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 103, 186802 (2009).


# Model-Free Learning of Optimal Beamformers for
Passive IRS-Assisted Sumrate Maximization

###### Abstract

Although Intelligent Reflective Surfaces (IRSs) are a cost-effective
technology promising high spectral efficiency in future wireless networks,
obtaining optimal IRS beamformers is a challenging problem with several
practical limitations. Assuming fully-passive, sensing-free IRS operation, we
introduce a new data-driven Zeroth-order Stochastic Gradient Ascent (ZoSGA)
algorithm for sumrate optimization in an IRS-aided downlink setting. ZoSGA
does not require access to channel model or network structure information, and
enables learning of optimal long-term IRS beamformers jointly with standard
short-term precoding, based only on conventional effective channel state
information. Supported by state-of-the-art (SOTA) convergence analysis,
detailed simulations confirm that ZoSGA exhibits SOTA empirical behavior as
well, consistently outperforming standard fully model-based baselines, in a
variety of scenarios.

Index Terms— Intelligent Reflecting Surfaces, Sumrate Maximization, Zeroth-
order Optimization, Model-Free Learning.

##  1 Introduction

Intelligent Reflective (or Reconfigurable Intelligent) Surfaces (IRSs or RISs)
are planar surfaces comprised of passive reflective elements with tunable
phase-shifts and amplitudes. Seen as beamforming arrays, IRSs present a
promising solution to mitigate sharp drops in Quality-of-Service (QoS) in the
absence of line of sight for highly directional mmWave signals, and beyond.
Such a typical scenario is visualized in Fig. 1, where users linked with an
Access Point (AP) are enjoying improved QoS due to the IRSs assisting the
network.

In IRS-aided communications, the goal is to optimally tune the IRS elements
along with other possible resources (such as AP precoders) to optimize a
certain system utility. A standard case is that of the weighted sumrate
utility in a MISO downlink scenario (see Fig. 1), where the goal is to
maximize the total downlink rate of a number of users (i.e., terminals)
actively serviced by an AP, while passively aided by one or multiple IRSs.
Then, the objective is to jointly optimize IRS amplitudes/phase-shifts and AP
precoders, subject to certain power constraints. While AP precoders are
usually continuous-valued, IRS phase-shifts can be either quantized [1] or
continuously varying [2, 3]. In this paper, we focus on the latter, i.e., on
optimizing for continuous IRS phase-shift variations, however strictly
assuming fully-passive IRS elements with no sensing capabilities or extra
hardware or scheduling requirements.

Recently, various deep learning driven methods have been proposed for IRS-
aided beamforming optimization. Offline learning approaches have been proposed
primarily to estimate the Channel State Information (CSI) using labeled
datasets with Function Approximators (FAs) [4, 5, 6, 7]. In contrast, Deep
Reinforcement Learning (DRL) methods have been used for joint beamforming
optimization, either with quantized [8, 9, 10] or continuous [11, 12] actions.
Though such end-to-end DRL methods do not require intermediate CSI estimates
or active IRS sensing, they do require FAs to at least approximate value
functions, jointly modeling the whole optimization task, including beamforming
controllers. On the one hand, choosing such FAs without explicit domain
knowledge usually results in increased problem complexity, non-
interpretability and lack of robustness, whereas incorporating specific domain
knowledge often results in overfitting, hindering versatility and
transferability to distinct environments. Further, as such approaches
primarily consider reactive IRS tuning (i.e., phase shifts depend on observed
CSI), they naturally incur increased power consumption for perpetual IRS
control.

Fig. 1: Realization/Concept of an IRS-aided Wireless Network.

To improve complexity, efficiency and robustness, model-based two-timescale
approaches have also been proposed, in which reactive (AP) precoding is
combined with static (i.e., non-reactive) IRS beamforming [13, 14, 15, 16,
17]. Although potentially effective, such methods rely on optimizing problem
surrogates and leverage detailed knowledge of channel models and spatial
network configurations, often requiring active sensing at the IRSs. Thus, they
are not applicable to different settings without complete remodeling of the
environment.

In this paper, we develop a new Zeroth-order Stochastic Gradient Ascent
(ZoSGA) two-timescale algorithm for weighted sumrate optimization in the MISO
downlink setting through purely data-driven learning of optimal IRS
beamformers, and joint reactive AP precoding via the standard WMMSE algorithm
[18]. ZoSGA follows the two-timescale paradigm while being completely FA- and
model-free, and relies on minimal zeroth-order reinforcement (i.e., system
probing), as well as user-experienced effective CSI (see Section 2),
conventionally available in mutli-user downlink beamfoming, regardless of the
number or spatial configuration of the IRSs. In fact, ZoSGA allows treating
the IRSs assisting the network as fully passive (though tunable) reflecting
elements, and applies to different networking scenarios seamlessly.

We analyze ZoSGA under mild regularity conditions and establish a state-of-
the-art (SOTA) convergence rate of order $\mathcal{O}(\sqrt{S}\epsilon^{-4})$,
where $S$ is the number of IRS tunable parameters and $\epsilon$ is the
suboptimality target. Then, we numerically demonstrate that ZoSGA also
achieves SOTA empirical performance in a variety of scenarios, substantially
outperforming the two-timescale method of [14], which is a standard model-
based baseline for the problem under study, assuming full knowledge of the
underlying channel model and spatial network configuration. Overall, ZoSGA
sets a new SOTA for beamforming optimization in IRS-aided communications.

Note: Proofs are omitted due to lack of space and will be included in a
journal paper currently under preparation. Fully reproducible code is
available at github.com/hassaanhashmi/zosga-irs.

##  2 Problem Formulation

Capitalizing on the standardized IRS-Aided communication setting depicted in
Fig. 1, our goal is to maximize the total downlink rate of $K$ users actively
serviced by an AP with $M$ antennas, while passively aided by one or multiple
IRSs, arbitrarily spatially placed. As explained in Section 1, we assume
dynamic (i.e., reactive) AP beamformers (also called precoders), while the IRS
beamformers are realistically viewed as static (i.e., non-reactive) elements,
tunable by a parameter vector $\bm{\theta}\in\mathbb{R}^{S}$, encompassing
whatever propagation feature is learnable on the IRSs present in the network.
For instance, $\bm{\theta}$ might refer to amplitudes, phases, or both, or
even realistic physical elements of an IRS, such as tunable varactor
capacitance elements, see, e.g., [2, 3]. We strictly make no sensing
assumptions on the IRSs, i.e., the IRSs are completely passive network
elements.

Each user $k=1,\ldots,K$ experiences a random effective channel denoted by
$\mathbf{h}_{k}\left(\bm{\theta},\omega\right)$, indexed by the IRS parameter
vector $\bm{\theta}$ as well as a state of nature $\omega\in\Omega$ describing
unobservable random propagation patterns for each value of $\bm{\theta}$. In
other words, $\mathbf{h}_{k}\left(\bm{\theta},\omega\right)$ is a random
channel field with spatial variable $\bm{\theta}$. We make the standard
assumption that the effective channels
$\mathbf{h}_{k}\left(\bm{\theta},\omega\right),k=1,\ldots,K$ are known to the
AP at the time of transmission [13, 14]. Note that the implementation
complexity of estimating effective channels in our setting is exactly the same
as that in conventional multi-user downlink beamforming (i.e., involving no
IRSs), regardless of the number and/or spatial configuration of the IRSs; no
extra hardware or customized scheduling schemes (such as those in [14]) are
required on either the AP or the IRSs assisting the network.

The QoS of user $k$ is measured by the corresponding SINR, i.e.,

$\displaystyle\begin{aligned}
\text{SINR}_{k}(\mathbf{W},\mathbf{h}_{k}(\bm{\theta},\omega))\triangleq\frac{\left|\mathbf{h}_{k}^{\mathsf{H}}(\bm{\theta},\omega)\mathbf{w}_{k}\right|^{2}}{{\sum}_{{j\in\mathbb{N}_{K}^{+}\setminus
k}}\left|\mathbf{h}^{\mathsf{H}}_{k}(\bm{\theta},\omega)\mathbf{w}_{j}\right|^{2}+\sigma^{2}_{k}},\end{aligned}$

where
$\mathbf{W}{=}\begin{bmatrix}\mathbf{w}_{1}&\mathbf{w}_{2}&\cdots&\mathbf{w}_{K}\end{bmatrix}\in\mathbb{C}^{M\times
K}$, $\mathbf{w}_{k}$ is a transmit precoding vector and $\sigma^{2}_{k}$ is
the noise variance for user $k$, respectively. Then, the weighted sumrate
utility of the network is defined as

${F}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))\triangleq\sum_{k=1}^{K}\alpha_{k}\log_{2}\left(1+\textnormal{SINR}_{k}\left(\mathbf{w},\mathbf{h}_{k}\left(\bm{\theta},\omega\right)\right)\right),$

with
$\mathbf{H}=\begin{bmatrix}\mathbf{h}_{1}\ldots\mathbf{h}_{K}\end{bmatrix}\in\mathbb{C}^{M\times
K}$, and $\alpha_{k}\geq 0$ the weight associated with user $k$. We are
interested in maximizing the sumrate of the network jointly by selecting
instantaneous-optimal dynamic AP precoders $\mathbf{W}$, and on-average-
optimal static IRS beamformers $\bm{\theta}$ [13, 14], i.e., we are interested
in the problem

$\boxed{\max_{\bm{\theta}\in\mathcal{K}}\mathbb{E}\bigg{\\{}\max_{\mathbf{W}:\|\mathbf{W}\|_{F}^{2}\leq
P}{F}\left(\mathbf{W},\mathbf{H}(\bm{\theta},\omega)\right)\bigg{\\}},}$ (1)

where $\|\cdot\|_{F}^{2}$ denotes the Frobenius norm, $P>0$ is a total power
budget at the AP, and $\mathcal{K}$ is a compact feasible set. Problem (1) is
an instance of a $2$-stage stochastic program [19]. Accordingly, we may
interpret optimization over $\bm{\theta}$ as the first-stage problem, where
the decision-maker optimizes the IRS parameters on-average and before the
actual effective channels are revealed, whereas optimization over $\mathbf{W}$
constitutes the second stage problem which is solved after randomness is
revealed to the decision-maker, therefore interpreting optimal AP beamforming
as (optimal) recourse actions.

$\displaystyle\begin{aligned}
\bm{D}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))=\sum_{k=1}^{K}\alpha_{k}\frac{\mathbf{z}^{\mathsf{H}}\left[\left(\sum_{j\in\mathbb{N}_{K}^{+}\setminus
k}\left|\mathbf{z}^{\mathsf{H}}\mathbf{w}_{j}\right|^{2}+\sigma^{2}_{k}\right)\mathbf{w}_{k}{\mathbf{w}_{k}}^{\mathsf{H}}-\left|\mathbf{z}^{\mathsf{H}}\mathbf{w}_{k}\right|^{2}\sum_{j\in\mathbb{N}_{K}^{+}\setminus
k}\mathbf{w}_{j}{\mathbf{w}_{j}}^{\mathsf{H}}\right]}{\ln(2)(\sum_{j\in\mathbb{N}_{K}^{+}\setminus
k}|\mathbf{z}^{\mathsf{H}}\mathbf{w}_{j}|^{2}+\sigma^{2}_{k})^{2}+|\mathbf{z}^{\mathsf{H}}\mathbf{w}_{k}|^{2}(\sum_{j\in\mathbb{N}_{K}^{+}\setminus
k}|\mathbf{z}^{\mathsf{H}}\mathbf{w}_{j}|^{2}+\sigma^{2}_{k})}\Bigg{|}_{\mathbf{z}=\mathbf{H}(\bm{\theta},\omega)}.\end{aligned}$
(D) $\displaystyle\begin{aligned}
\dfrac{\mathbf{G}_{\mu}(\bm{\theta},\omega,\mathbf{u})}{2}\triangleq&\left(\frac{\Re\left(\mathbf{\Delta}_{\mu}(\bm{\theta},\omega,\mathbf{u})\right)}{2\mu}\mathbf{u}^{\mathsf{T}}\right)^{\mathsf{T}}\left(\Re\left(\bm{D}(\mathbf{W}^{*},\mathbf{H}(\bm{\theta},\omega))\right)\right)^{\mathsf{T}}+\left(\frac{\Im\left(\mathbf{\Delta}_{\mu}(\bm{\theta},\omega,\mathbf{u})\right)}{2\mu}\mathbf{u}^{\mathsf{T}}\right)^{\mathsf{T}}\left(\Re\left(j\bm{D}(\mathbf{W}^{*},\mathbf{H}(\bm{\theta},\omega))\right)\right)^{\mathsf{T}}.\end{aligned}$
(G)

##  3 An inner-outer optimization method

We would like to devise a gradient-inspired method for tackling problem (1),
which circumvents the need for accessing first-order gradient information of
the effective channel function $\mathbf{H}(\cdot,\omega)$; in principle, this
is unknown. To that end, we derive an inner-outer scheme for the solution of
the two-stage problem in (1). At each (outer) iteration (i.e., for each state
of nature $\omega$), we employ an oracle to find an approximate solution to
the inner problem, which is then utilized to derive a model-free gradient
approximation for the outer problem using zeroth-order reinforcement (i.e.,
probing) on the involved effective channel. This approach completely bypasses
the need of a model for the effective channel, and offers great flexibility as
to how the inner problem can be approximately solved.

###  3.1 Tackling the Inner Problem (AP Precoding)

To obtain an approximate solution for the (nonconvex) inner maximization
problem, we heuristically employ the well-known weighted minimum mean squared
error (WMMSE) algorithm [18], simultaneously implementing AP precoding. For
technical purposes, we also introduce the following assumption.

######  Assumption 1.

Given any $\bm{\theta}\in\mathcal{K}$ and almost every (a.e.)
$\omega\in\Omega$, we have access to an oracle that yields an optimal solution
$\bm{W}^{*}(\bm{\theta},\omega)\in\arg\max_{\mathbf{W}\colon\|\mathbf{W}\|_{F}^{2}\leq
P}{F}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))$.

In practice Assumption 1 may be restrictive but, as verified in Section 5, the
empirical effectiveness of the proposed algorithm is not hindered.
Nevertheless, this assumption is crucial for the derivation of our algorithm
and its subsequent theoretical analysis.

###  3.2 Tackling the Outer Problem (IRS Beamforming)

Assuming the availability of an oracle providing an optimal solution to the
inner problem at any $(\bm{\theta},\omega)$, say
$\mathbf{W}^{*}(\bm{\theta},\omega)$ (in practice obtained approximately via
WMMSE), we can write problem (1) as

$\begin{split}\max_{\bm{\theta}\in\mathcal{K}}&\
\mathbb{E}\left\\{{F}\left(\mathbf{W}^{*}(\bm{\theta},\omega),\mathbf{H}(\bm{\theta},\omega)\right)\right\\}.\end{split}$
(2)

As we would like to derive a gradient-ascent-like scheme for solving (2), let
us make the following assumptions on the effective channel.

######  Assumption 2.

For almost all $\omega$, the channel function $\mathbf{H}(\cdot,\omega)$ is
uniformly bounded, twice continuously-differentiable,
$L_{\mathbf{h},0}$-Lipschitz, with $L_{\mathbf{h},1}$-Lipschitz gradients.

We note that Assumption 2 imposes regularity conditions mainly required for
the grounded development of our optimization scheme and for its convergence
analysis (later in Section 4). Also, observe that the boundedness assumption
is natural, since (IRS-aided) wireless channels are always bounded in
practice. While we usually have no information on the analytical properties of
the effective channel, Assumption 2 is easily satisfied in widely used channel
models of IRS-aided systems, see, e.g., [14] as well as Section 5.

From Assumptions 1–2, and the compactness of $\mathcal{K}$, we obtain

$\begin{split}&\nabla_{\bm{\theta}}\mathbb{E}\left\\{{F}\left(\mathbf{W}^{*}(\bm{\theta},\omega),\mathbf{H}(\bm{\theta},\omega)\right)\right\\}\\\
&\quad=\mathbb{E}\left\\{\nabla_{\bm{\theta}}{F}\left(\mathbf{W}^{*}(\bm{\theta},\omega),\mathbf{H}(\bm{\theta},\omega)\right)\right\\}\\\
&\quad=\mathbb{E}\left\\{\nabla_{\bm{\theta}}{F}\left(\mathbf{W},\mathbf{H}\left(\bm{\theta},\omega\right)\right)\mid_{\mathbf{W}=\mathbf{W}^{*}(\bm{\theta},\omega)}\right\\},\end{split}$
(3)

where we used [19, Theorem 7.44], and the _implicit function theorem_ (e.g.,
see [20]). Next, we derive the gradient of
${F}(\mathbf{W},\mathbf{H}(\cdot,\omega))$.

####  3.2.1 Gradient Representation for the Outer Problem

We notice that ${F}(\mathbf{W},\cdot)$ takes a complex input and hence we need
an appropriate generalization of the gradient for it. To that end, we utilize
the so-called _Wirtinger calculus_ (see [21]). In the following lemma we
derive the full compositional gradient of
${F}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))$ by following the developments
in [21, Section 4].

######  Lemma 1.

For every $\bm{\theta}\in\mathcal{K}$, $\mathbf{W}\in\mathbb{C}^{M\times K}$
and $\omega\in\Omega$, the gradient of
${F}\left(\mathbf{W},\mathbf{H}(\bm{\theta},\omega)\right)$ with respect to
$\bm{\theta}$ reads

$\begin{split}\nabla_{\bm{\theta}}&{F}\left(\mathbf{W},\mathbf{H}(\bm{\theta},\omega)\right)\\\
&\quad=2\nabla_{\bm{\theta}}\Re\left(\mathbf{H}(\bm{\theta},\omega)\right)\left(\Re\left(\bm{D}\left(\mathbf{W},\mathbf{H}(\bm{\theta},\omega)\right)\right)\right)^{\mathsf{T}}\\\
&\qquad+2\nabla_{\bm{\theta}}\Im\left(\mathbf{H}(\bm{\theta},\omega)\right)\left(\Re\left(j\bm{D}\left(\mathbf{W},\mathbf{H}(\bm{\theta},\omega)\right)\right)\right)^{\mathsf{T}},\end{split}$

where $\bm{D}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))$ is given in (D).
Additionally, there exists a constant $B_{F}>0$ such that
$\left\|\bm{D}(\mathbf{W},\mathbf{H}(\bm{\theta},\omega))\right\|\leq B_{F}$
for every $(\mathbf{W},\bm{\theta})\in\mathbb{C}^{M\times K}\times\mathcal{K}$
and a.e. $\omega\in\Omega$, where $\mathbf{W}$ is feasible for (1).

####  3.2.2 Model-Free Jacobian Approximation of the Effective Channel

Assuming the lack of availability of any first-order information of
$\mathbf{H}(\cdot,\omega)$ for any $\omega\in\Omega$, we will employ a zeroth-
order scheme in order to obtain a Jacobian estimate of
$\nabla_{\bm{\theta}}\mathbf{H}(\bm{\theta},\omega)$, using which we can solve
(2) via a stochastic gradient ascent scheme. The proposed method will be based
on Jacobian estimates arising from a two-point stochastic evaluation of
$\mathbf{H}(\cdot,\omega)$ (similar to, among others, [22, 23, 24]). From
Assumption 2, we can first write

$\nabla_{\bm{\theta}}\mathbf{H}(\bm{\theta},\omega)=\nabla_{\bm{\theta}}\Re\left(\mathbf{H}(\bm{\theta},\omega)\right)+j\nabla_{\bm{\theta}}\Im\left(\mathbf{H}(\bm{\theta},\omega)\right).$

We approximate this gradient using only function evaluations of
$\mathbf{H}(\cdot,\omega)$ (i.e. by minimally probing the network). For each
$\omega\in\Omega$, let
$\mathbf{u}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}\right)$. Given a
smoothing parameter $\mu>0$, we define

$\begin{split}\nabla^{\mu}_{\bm{\theta}}\mathbf{H}(\bm{\theta},\omega)\triangleq&\
\mathbb{E}\left\\{\frac{\mathbf{H}\left(\bm{\theta}+\mu\mathbf{u},\omega\right)-\mathbf{H}\left(\bm{\theta}-\mu\mathbf{u},\omega\right)}{2\mu}\mathbf{u}^{\mathsf{T}}\right\\}^{\mathsf{T}}.\end{split}$
(4)

Using Lemma 1, we substitute the Jacobian of the effective channel via the
zeroth-order approximation given in (4) to obtain the zeroth-order
compositional quasi-gradient (an approximation of (2)) as

$\begin{split}{\hat{\nabla}}_{\bm{\theta}}^{\mu}&{F}\left(\mathbf{W}^{*},\mathbf{H}(\bm{\theta},\omega)\right)\\\
&\quad=2\nabla^{\mu}_{\bm{\theta}}\Re\left(\mathbf{H}(\bm{\theta},\omega)\right)\left(\Re\left(\bm{D}\left(\mathbf{W}^{*},\mathbf{H}(\bm{\theta},\omega)\right)\right)\right)^{\mathsf{T}}\\\
&\qquad+2\nabla^{\mu}_{\bm{\theta}}\Im\left(\mathbf{H}(\bm{\theta},\omega)\right)\left(\Re\left(j\bm{D}\left(\mathbf{W}^{*},\mathbf{H}(\bm{\theta},\omega)\right)\right)\right)^{\mathsf{T}}.\end{split}$

We are now able to derive the proposed inner-outer scheme for the solution of
problem (1). At every iteration $t$, we draw independent and identically
distributed (i.i.d.) samples $\omega^{t+1}$ and $\mathbf{u}^{t+1}$, and
compute an “optimal" inner solution
$\mathbf{W}^{*}(\bm{\theta}^{t},\omega^{t+1})$, as described in Section 3.1.
Then, by probing the system twice, we can compute a zeroth-order stochastic
gradient estimator as shown in (G), where
$\mathbf{\Delta}_{\mu}(\bm{\theta},\omega,\mathbf{u})\triangleq\mathbf{H}(\bm{\theta}+\mu\mathbf{u},\omega)-\mathbf{H}(\bm{\theta}-\mu\mathbf{u},\omega)$.
The proposed Zeroth-order Stochastic Gradient Ascent (ZoSGA) method with WMMSE
is provided in Algorithm 1.

1 Initialize $\mathbf{W}^{0}$, $\bm{\theta}^{0}$, $\eta^{0}$, $\mu$.

2 for  _$t=0,\ldots,T$_ do

3 Sample (i.i.d.) $(\mathbf{u}^{t+1},\omega^{t+1})$.

4 Find $\mathbf{W}^{*}(\bm{\theta}^{t},\omega^{t+1})$ using WMMSE [18].

5 Compute $\mathbf{G}_{\mu}(\bm{\theta}^{t},\omega^{t+1},\mathbf{u}^{t+1})$
using (G) (by probing).

6 Choose $\eta^{t+1}$ and update
$\bm{\theta}^{t+1}=\Pi_{\mathcal{K}}\left(\bm{\theta}^{t}+\eta^{t+1}\mathbf{G}_{\mu}(\bm{\theta}^{t},\omega^{t+1},\mathbf{u}^{t+1})\right),$
where $\Pi_{\mathcal{K}}(\cdot)$ is the Euclidean projection onto
$\mathcal{K}$.

7 end for

Return $\bm{\theta}^{t^{*}}$, where
$\mathbb{P}(t^{*}=t)=\frac{\eta^{n}}{\sum_{i=0}^{T}\eta^{i}}$.

Algorithm 1 ZoSGA with WMMSE

(a)

(b)

(c)

Fig. 2: Average sumrates achieved by WMMSE [18] (with random IRS phase
shifts), TTS-SSCO [14], and the proposed ZoSGA algorithm, with only IRS 1
present: (a) with fixed parameters $\beta_{*}$ and $r_{*}$, and (b) relative
to Rician factors $\beta_{AI}=\beta_{Iu}$; (c) Average sumrates achieved by
WMMSE and ZoSGA with IRS 1+2 present and fixed parameters $\beta_{*}$ and
$r_{*}$. (AA: Adjustable Amplitude; UA: Unit Amplitude)

##  4 Convergence analysis

In this section we provide a brief account of the convergence properties of
Algortithm 1. For convenience let
$\varphi(\cdot)\triangleq-\mathbb{E}\\{F_{\omega}^{*}(\cdot)\\}+\delta_{\mathcal{K}}(\cdot)$,
where $F_{\omega}^{*}(\cdot)\triangleq
F\left(\mathbf{W}^{*}(\cdot,\omega),\mathbf{H}(\cdot,\omega)\right)$ is the
(sample) objective of problem (2) and
$\delta_{\mathcal{K}}(\cdot)\in\\{0,\infty\\}$ is the indicator of the set
$\mathcal{K}$. The following lemma is crucial for the analysis of Algorithm 1.

######  Lemma 2.

Let Assumptions 1, 2 be in effect, and let
$\widehat{\mathcal{K}}\supset\mathcal{K}$ be a convex compact set. Then, for
a.e. $\omega$ there exists $\rho(\omega)>0$ such that $-F_{\omega}^{*}$ is
$\rho(\omega)$-weakly convex on $\widehat{\mathcal{K}}$, i.e.,
$-F_{\omega}^{*}(\cdot)+(\rho(\omega)/2)\|\cdot\|^{2}$ is convex on
$\widehat{\mathcal{K}}$.

Given some $\lambda>0$, we define the _Moreau envelope_ of $\varphi$ as

$\varphi^{\lambda}(\bm{y})\triangleq\min_{\bm{\theta}}\left\\{\varphi(\bm{\theta})+\frac{1}{2\lambda}\|\bm{y}-\bm{\theta}\|^{2}\right\\}.$

A near-stationary point for the Moreau envelope $\varphi^{\lambda}$ is close
to a near-stationary solution of problem (1), if
$-\mathbb{E}\\{F_{\omega}^{*}(\cdot)\\}$ is weakly convex (see [25, Section
2.2] and note, from Lemma 2, that $-\mathbb{E}\\{F_{\omega}^{*}(\cdot)\\}$ is
weakly convex if $\rho(\omega)$ is integrable). In what follows, we provide a
convergence rate of Algorithm 1 in terms of the gradient of the Moreau
envelope.

######  Theorem 3.

Let Assumptions 1, 2 be in effect and assume that $\mathcal{K}$ is convex and
compact, $\rho=\mathbb{E}\\{\rho(\omega)\\}<\infty$, while the process
$\\{\omega^{t}\\}_{t\in\mathbb{N}}$ is i.i.d.. Let
$\\{\bm{\theta}^{t}\\}_{t=0}^{T}$ be the sequence of iterates generated by
Algorithm 1, $\Delta_{\mathcal{K}}>0$ be the diameter of $\mathcal{K}$, and
set

$\eta^{t}=\sqrt{\frac{\Delta_{\varphi}}{4\rho
B_{F}^{2}L_{\mathbf{h},0}^{2}(S^{2}+2S)(T+1)}},\qquad\textnormal{for all
}t\geq 0,$

for some
$\Delta_{\varphi}\geq\varphi^{1/(2\rho)}(\bm{\theta}^{0})-\min_{\bm{\theta}}\varphi(\bm{\theta})$.
Then, it is true that

$\begin{split}&\mathbb{E}\big{\\{}\big{\|}\nabla\varphi^{1/(2\rho)}\big{(}\bm{\theta}^{t^{*}}\big{)}\big{\|}^{2}_{2}\big{\\}}\\\
&\leq\ 8\Bigg{(}\sqrt{\frac{\Delta_{\varphi}\rho
B_{F}^{2}L_{\mathbf{h},0}^{2}(S^{2}+2S)}{T+1}}+\mu\rho\Delta_{\mathcal{K}}B_{F}L_{\mathbf{h},1}\sqrt{SMK}\Bigg{)}.\end{split}$

######  Remark 1.

Note that choosing $\mu=\mathcal{O}(1/\sqrt{(MKT)})$ yields that
$\mathbb{E}\big{\\{}\big{\|}\nabla\varphi^{1/(2\rho)}\big{(}\bm{\theta}^{t^{*}}\big{)}\big{\|}_{2}\big{\\}}\leq\epsilon$,
after $\mathcal{O}(\sqrt{S}\epsilon^{-4})$ iterations.

##  5 Simulations

We consider an IRS-aided MISO downlink wireless network with Rician fading
channels $\mathbf{G}_{i}\in\mathbb{C}^{N{\times}M}$,
$\mathbf{h}^{i}_{r,k}\in\mathbb{C}^{N}$ and
$\mathbf{h}_{d,k}\in\mathbb{C}^{M}$, respectively, as shown in Fig. 3, where
$i\in\\{1,2\\}$. Following [14], these channels are defined for each user $k$
as

$\displaystyle\begin{aligned}
\mathbf{h}_{r,k}^{i}&{\triangleq}\sqrt{\beta_{Iu}{/}(1{+}\beta_{Iu})}\check{\mathbf{v}}^{i}_{r,k}{+}\sqrt{1{/}(1{+}\beta_{Iu})}\mathbf{\Phi}_{r,k}^{1/2}\mathbf{v}_{r,k}^{i},\\\
\mathbf{G}_{i}&{\triangleq}\sqrt{\beta_{AI}{/}(1{+}\beta_{AI})}\check{\mathbf{F}}^{i}{+}\sqrt{1{/}(1{+}\beta_{AI})}\mathbf{\Phi}_{r}^{1/2}\mathbf{F}^{i}\mathbf{\Phi}_{d}^{1/2},\\\
\mathbf{h}_{d,k}&{\triangleq}\sqrt{\beta_{Au}{/}(1{+}\beta_{Au})}\check{\mathbf{v}}_{d,k}{+}\sqrt{1{/}(1{+}\beta_{Au})}\mathbf{\Phi}_{d}^{1/2}\mathbf{v}_{d,k},\\\
\end{aligned}$ (5)

where the entries of
$\mathbf{H}_{\text{S-CSI}}{\triangleq}\\{\check{\mathbf{v}}^{i}_{r,k}{,}\check{\mathbf{F}}^{i}{,}\check{\mathbf{v}}_{d,k}\\}$
and $\mathbf{H}_{\text{I-CSI}}{\triangleq}$
$\\{\mathbf{v}^{i}_{r,k}{,}\mathbf{F}^{i}{,}$ $\mathbf{v}_{d,k}\\}$
($i\in\\{1,2\\}$) are all sampled independently from
$\mathcal{CN}(\bf{0},\bf{I})$, with the former being sampled only once per
simulation, where a simulation is defined as an ensemble containing a number
of sequential i.i.d. realizations of "instantaneous" CSI
$\mathbf{H}_{\text{I-CSI}}$ (given a fixed value for "statistical" CSI
$\mathbf{H}_{\text{S-CSI}}$ for each simulation).

Each IRS has $N{=}N_{h}{\times}N_{v}$ phase-shift elements (see Fig. 3). Thus,
$\bm{\theta}_{i}{=}\begin{bmatrix}\bm{\phi}_{i}^{\mathsf{T}}&\mathbf{A}_{i}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
where $\bm{\phi}_{i}\in[-\pi,\pi]^{N}$ and $\mathbf{A}_{i}\in[0,1]^{N}$ are
phases and amplitudes of the IRS elements, respectively, with $i\in\\{1,2\\}$.
The scalars $\beta_{Iu}$, $\beta_{AI}$, and $\beta_{Au}$, are link-specific
Rician factors, while $\mathbf{\Phi}_{r,k}\in\mathbb{C}^{N{\times}N}$,
$\mathbf{\Phi}_{r}\in\mathbb{C}^{N{\times}N}$ and
$\mathbf{\Phi}_{d}\in\mathbb{C}^{M{\times}M}$ are spatial correlation matrices
following the exponentially decaying correlation model of [14, 15], with
corresponding correlation coefficients $r_{r,k}\in(0,1)$, $r_{r}\in(0,1)$ and
$r_{d}\in(0,1)$. Lastly, the distance-dependent path-loss for each link is
modeled as $L=\sqrt{C_{0}{d}^{-\alpha}}$, where $d$ is distance (in meters),
$\alpha$ is the corresponding path loss exponent, and $C_{0}{=}$-30dBm.

Throughout the simulations we consider a power allocation of $P{=}$5dBm, a
noise variance $\sigma^{2}{=}$ -80dBm, and run WMMSE for 20 iterations to
optimize the precoding vectors $\mathbf{W}$ for each instance of the received
effective channels $\mathbf{h}_{k}(\bm{\theta},\omega)$, each of which is
given as

$\displaystyle\mathbf{h}_{k}(\bm{\theta},\mathbf{\omega})=\sum_{i=1}^{2}\underbrace{\mathbf{G}^{\mathsf{H}}_{i}\text{diag}(\mathbf{A}_{i}\circ
e^{-j{\bm{\phi}}_{i}})\mathbf{h}^{i}_{r,k}}_{\bm{\theta}_{i}\text{-reflected
link}}+\underbrace{\mathbf{h}_{d,k}}_{\text{LoS link}},$ (6)

for all $k\in\mathbb{N}^{+}_{4}$, where
$\bm{\theta}{=}(\bm{\theta}_{1},\bm{\theta}_{2})$,
$\omega=\\{\mathbf{G}_{1},\mathbf{G}_{2},\mathbf{h}^{1}_{r,k},\mathbf{h}^{2}_{r,k},\mathbf{h}_{d,k},$
$k=1,\ldots,K\\}$ and, for simplicity, we ignore IRS-to-IRS links. All
presented results are averaged over 2000 unique simulations.

Fig. 3: Experimental IRS-Aided Network Configuration.

First, we compare the proposed ZoSGA against TTS-SSCO [14] for a network with
one IRS as TTS-SSCO requires an exact model to optimize the IRS phase-shift
elements (see [14, eq. (35)]). We let $\mu{=}10^{-12}$, and choose separate
step-sizes $\eta^{0}_{\phi}=0.4$ and $\eta^{0}_{A}=0.01$ for $\bm{\phi}$ and
$\mathbf{A}$, respectively scaled by ${0.9972}^{t}$ for
$t\in\mathbb{N}^{+}_{10^{3}}$, keeping them constant afterwards. The
parameters for TTS-SSCO are provided in [14, Section V]. From Fig. 2(a), we
observe that ZoSGA substantially outperforms TTS-SSCO solely on the basis of
effective CSI, while having no access to the statistical model of the channel
or the spatial configuration of the system; this is in sharp contrast to TTS-
SSCO. While TTS-SSCO converges faster, it does so by internally sampling I-CSI
and finding corresponding optimal precoding vectors (via WMMSE) 10 times per
iteration. In Fig. 2(b), we see that, for spatially uncorrelated channels, the
relative gain of ZoSGA in the achievable sumrate increases with respect to the
Rician factors pertaining to $\bm{\theta}_{1}$-reflected links, i.e., as we
move from I-CSI to S-CSI dominated channels. We also observe that WMMSE with a
randomized IRS is insensitive to changes in Rician factor; this is expected.

Lastly, we use both IRSs and optimize their phase-shift elements using
identical stepsizes for ZoSGA as stated above. Our results in Fig. 2(c) show
that ZoSGA not only scales well to unknown system/channel models, but that it
is also robust to the choices of $\eta_{\phi}$ and $\eta_{A}$. The results
also demonstrate the better performance gains when optimizing
$\bm{\theta}_{1}$ versus optimizing the more distant $\bm{\theta}_{2}$.

##  6 Conclusion

We introduced ZoSGA, a truly model-free method, to optimize fully-passive IRS
phase-shift elements for a given QoS metric. We proved state-of-the-art
convergence rates for ZoSGA under standard assumptions, and empirically
demonstrated its efficiency on a benchmark IRS-aided network. In fact, ZoSGA
outperforms the current SOTA by learning (near-)optimal passive IRS
beamformers solely on the basis of conventional effective CSI, in the absence
of channel models and spatial network configuration information.

## References

  * [1] H. Kamoda, T. Iwasaki, J. Tsumochi, T. Kuki, and O. Hashimoto,  “60-ghz electronically reconfigurable large reflectarray using single-bit phase shifters,”  IEEE transactions on antennas and propagation, vol. 59, no. 7, pp. 2524–2531, 2011.
  * [2] J. Zhao, Q. Cheng, J. Chen, M. Q. Qi, W. X. Jiang, and T. J. Cui,  “A tunable metamaterial absorber using varactor diodes,”  New Journal of Physics, vol. 15, no. 4, pp. 043049, 2013.
  * [3] A. Araghi, M. Khalily, M. Safaei, A. Bagheri, V. Singh, F. Wang, and R. Tafazolli,  “Reconfigurable intelligent surface (ris) in the sub-6 ghz band: Design, implementation, and real-world demonstration,”  IEEE Access, vol. 10, pp. 2646–2655, 2022.
  * [4] A. Taha, M. Alrabeiah, and A. Alkhateeb,  “Deep learning for large intelligent surfaces in millimeter wave and massive mimo systems,”  in 2019 IEEE Global communications conference (GLOBECOM). IEEE, 2019, pp. 1–6.
  * [5] E. Balevi, A. Doshi, A. Jalal, A. Dimakis, and J. G. Andrews,  “High dimensional channel estimation using deep generative networks,”  IEEE Journal on Selected Areas in Communications, vol. 39, no. 1, pp. 18–30, 2020.
  * [6] B. Yang, X. Cao, C. Huang, C. Yuen, L. Qian, and M. Di Renzo,  “Intelligent spectrum learning for wireless networks with reconfigurable intelligent surfaces,”  IEEE Transactions on Vehicular Technology, vol. 70, no. 4, pp. 3920–3925, 2021.
  * [7] S. Zhang, S. Zhang, F. Gao, J. Ma, and O. A. Dobre,  “Deep learning optimized sparse antenna activation for reconfigurable intelligent surface assisted communication,”  IEEE Transactions on Communications, vol. 69, no. 10, pp. 6691–6705, 2021.
  * [8] F. B. Mismar, B. L. Evans, and A. Alkhateeb,  “Deep reinforcement learning for 5g networks: Joint beamforming, power control, and interference coordination,”  IEEE Transactions on Communications, vol. 68, no. 3, pp. 1581–1592, 2019.
  * [9] A. Taha, Y. Zhang, F. B. Mismar, and A. Alkhateeb,  “Deep reinforcement learning for intelligent reflecting surfaces: Towards standalone operation,”  in 2020 IEEE 21st international workshop on signal processing advances in wireless communications (SPAWC). IEEE, 2020, pp. 1–5.
  * [10] H. Yang, Z. Xiong, J. Zhao, D. Niyato, L. Xiao, and Q. Wu,  “Deep reinforcement learning-based intelligent reflecting surface for secure wireless communications,”  IEEE Transactions on Wireless Communications, vol. 20, no. 1, pp. 375–388, 2020.
  * [11] C. Huang, R. Mo, and C. Yuen,  “Reconfigurable intelligent surface assisted multiuser miso systems exploiting deep reinforcement learning,”  IEEE Journal on Selected Areas in Communications, vol. 38, no. 8, pp. 1839–1850, 2020.
  * [12] Z. Yang, Y. Liu, Y. Chen, and J. T. Zhou,  “Deep reinforcement learning for ris-aided non-orthogonal multiple access downlink networks,”  in GLOBECOM 2020-2020 IEEE Global Communications Conference. IEEE, 2020, pp. 1–6.
  * [13] H. Guo, Y.-C. Liang, J. Chen, and E. G. Larsson,  “Weighted sum-rate maximization for reconfigurable intelligent surface aided wireless networks,”  IEEE Transactions on Wireless Communications, vol. 19, no. 5, pp. 3064–3076, 2020.
  * [14] M.-M. Zhao, Q. Wu, M.-J. Zhao, and R. Zhang,  “Intelligent reflecting surface enhanced wireless networks: Two-timescale beamforming optimization,”  IEEE Transactions on Wireless Communications, vol. 20, no. 1, pp. 2–17, 2020.
  * [15] M.-M. Zhao, A. Liu, Y. Wan, and R. Zhang,  “Two-timescale beamforming optimization for intelligent reflecting surface aided multiuser communication with qos constraints,”  IEEE Transactions on Wireless Communications, vol. 20, no. 9, pp. 6179–6194, 2021.
  * [16] Z. Yang, M. Chen, W. Saad, W. Xu, M. Shikh-Bahaei, H. V. Poor, and S. Cui,  “Energy-efficient wireless communications with distributed reconfigurable intelligent surfaces,”  IEEE Transactions on Wireless Communications, vol. 21, no. 1, pp. 665–679, 2021.
  * [17] Y. Liu, Q. Hu, Y. Cai, G. Yu, and G. Y. Li,  “Deep-unfolding beamforming for intelligent reflecting surface assisted full-duplex systems,”  IEEE Transactions on Wireless Communications, 2021.
  * [18] Q. Shi, M. Razaviyayn, Z. Q. Luo, and C. He,  “An iteratively weighted mmse approach to distributed sum-utility maximization for a mimo interfering broadcast channel,”  IEEE Transactions on Signal Processing, vol. 59, pp. 4331–4340, 9 2011.
  * [19] A. Shapiro, D. Dentcheva, and A. Ruszczyński,  Lectures on Stochastic Programming: Modeling and Theory,  MOS-SIAM Series on Optimization. SIAM & Mathematical Optimization Society, Philadelphia, 2014.
  * [20] A. L. Dontchev and R. T. Rockafellar,  Implicit Functions and Solution Mappings,  Springer Monographs in Mathematics. Springer New York, NY, 2009.
  * [21] K. Kreutz-Delgado,  “The complex gradient operator and the $\mathbb{C}\mathbb{R}$-calculus,”  arXiv preprint arXiv:0906.4835, 2009.
  * [22] D. S. Kalogerias and W. B. Powell,  “Zeroth-order stochastic compositional algorithms for risk-aware learning,”  SIAM Journal on Optimization, vol. 32, no. 2, 2022.
  * [23] Y. Nesterov and V. Spokoiny,  “Random gradient-free minimization of convex functions,”  Foundations of Computational Mathematics, vol. 17, pp. 527–566, 2017.
  * [24] S. Pougkakiotis and D. S. Kalogerias,  “A zeroth-order proximal stochastic gradient method for weakly convex stochastic optimization,”  arXiv preprint arXiv:2205.01633, 2022.
  * [25] D. Davis and D. Drusvyatskiy,  “Stocahstic model-based minimization of weakly convex functions,”  SIAM Journal on Optimization, vol. 29, no. 1, pp. 207–239, 2019\.