# On the Stability of Random Multiple Access with Feedback Exploitation and
Queue Priority

Karim G. Seddik
Electronics Engineering Department, American University in Cairo AUC Avenue,
New Cairo 11835, Egypt. email<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the stability of two interacting queues under random
multiple access in which the queues leverage the feedback information. We
derive the stability region under random multiple access where one of the two
queues exploits the feedback information and backs off under negative
acknowledgement (NACK) and the other, higher priority, queue will access the
channel with probability one. We characterize the stability region of this
feedback-based random access protocol and prove that this derived stability
region encloses the stability region of the conventional random access (RA)
scheme that does not exploit the feedback information.

##  I Introduction

The stability of interacting queues has been extensively considered in
literature. Several works have considered the characterization of the
stability region of interacting queues under random access protocols. The
stability region is characterized for the case $M=2$ and $M=3$ interacting
queues as well as the case of $M>3$ with symmetric arrivals. The stability
region for the general case of $M>3$ with asymmetric arrivals is still an open
problem and only inner achievable bounds are known.

Recently, many papers have considered the problem of interacting queues in
different contexts. For example, [1] considers the problem of interacting
queues in a TDMA system where a relay is used to help the source nodes in
forwarding their lost packets. In [2], the stability of interacting queues
under a random access protocol in the context of Cognitive Radio Network was
derived. In [3], the stability region of two interacting queues under random
access protocol where the two queues harvest energy was characterized. Other
works can be found in [4, 5], where derivations of the stability regions in
the context of different cognitive radio networks were considered.

In this paper, we derive the stability region of a two-queue random access
(RA) protocol with priorities. The queues will apply the conventional RA
protocol but in the case of packet loss due to collision the two queues will
exploit the feedback information to provide some level of coordination. We set
a priority to one of the two queues as follows. In the case of a negative
acknowledgement, the queue with the higher priority will attempt transmission
in the following time slot with probability one and the other queue will back
off to allow for collision-free transmission of the higher priority queue.
Clearly, this will enhance the service rate for the higher priority queue but
more interestingly it will also improve the service rate for the other, less
priority queue as will be explained later. We derive an expression for the
boundary of the stability region and prove that the RA with priority scheme
encloses the stability region of the conventional RA scheme.

To the best of our knowledge, the problem of characterizing the stability
region of the random access protocol with feedback leveraging has not been
considered before. We will characterize the stable arrival rates region and
prove that it contains that of the conventional random multiple access scheme
(with no feedback exploitation).

The rest of the paper is organized as follows. The system model is presented
in Section II. The performance of the proposed scheme is investigated in
Section III. The paper is concluded in Section IV. We have moved most of the
proofs to the appendices to preserve the flow of ideas in the paper.

##  II System Model

The system model is shown in Fig. 1. We consider the case of two interacting
packet queues, namely $Q_{1}$ and $Q_{2}$. $Q_{1}$ and $Q_{2}$ have infinite
buffers for storing fixed length packets. The channel is slotted in time and
any slot duration equals one packet transmission time. The arrival processes
at the two queues, $Q_{1}$ and $Q_{2}$, are modeled as Bernoulli arrival
processes with means $\lambda_{1}$ and $\lambda_{2}$, respectively [3]. Under
our system model assumptions, the average arrival rates are $\lambda_{1}$ and
$\lambda_{2}$ packets per time slot, and are bounded as $0<\lambda_{i}<1$,
$i=1,2$111The maximum service rate in our model is 1 packet/slot, since the
slot duration equals one packet transmission time, then the arrival rates must
be less than 1 otherwise the system will be unstable [3].. We can assume that
the packets arrive at the start of the time slot.

Figure 1: The system model.

The channel is modeled as a collision channel, where packet loss results only
in the case of simultaneous transmissions from the two queues. If only one
queue attempts transmitting at a given time slot, the packet is considered to
be correctly received [6, 3]. In the random access phase, the first queue
accesses the channel with probability $p_{1}$ whenever it has packets to send
and the second queue will access the channel with probability $p_{2}$ whenever
it has packets to send. If at any time slot some queue is empty, it will not
attempt any channel access.

In this paper, we will consider the use of the feedback information that is
leveraged at the queues in the case of collision. In the conventional random
multiple access system and in the case of collision, the collided packets stay
on the head of the queues and retransmissions are attempted employing the same
random multiple access scheme. In this paper, we consider a system where the
feedback information is leveraged at the queues and a priority is set to the
first queue; in the next time slot after collision, queue 2 ($Q_{2}$) will
back off and queue 1 ($Q_{1}$) will retransmit its collided packet to allow
for collision-free transmission of $Q_{1}$; after that the two queues return
to the conventional random multiple access scheme. The priority set to queue 1
can be due to some quality of service (QoS) requirement that is different from
the QoS requirement of queue 2. The interesting result is that although the
feedback will enhance the service of queue 1 by setting a higher priority to
it, the service will be enhanced as well for queue 2 as will be explained
later.

##  III The Stability Region for the Feedback-Based Random Access Protocol
with Priorities

In this section, we will characterize the stability region for the feedback-
based random access scheme. Stability can be loosely defined as having a
certain quantity of interest kept bounded. In our case, we are interested in
the queue size being bounded. For an irreducible and aperiodic Markov chain
with countable number of states, the chain is stable if and only if it is
positive recurrent, which implies the existence of its stationary
distribution. For a rigorous definition of stability under more general
scenarios see [6] and [7].

If the arrival and service processes of a queueing system are strictly
stationary, then one can apply Loynes’s theorem to check for stability
conditions [8]. This theorem states that if the arrival process and the
service process of a queueing system are strictly stationary, and the average
arrival rate is less than the average service rate, then the queue is stable,
otherwise it is unstable.

Characterizing the stability region will be a difficult problem due to the
interaction of the two queues and due to the fact that the service for one
queue will depend on the state of the other queue. We will consider the use of
the Dominant System concept that was proposed in [6] to characterize the
stability region of the conventional RA scheme. We will define two dominant
systems tailored to match our feedback-based random access scheme and in each
of the two systems we will determine the boundaries of the stability region.

###  III-A Dominant System 1

In any dominant system, we define a system that “stochastically dominates” our
system, that is the queues lengths in the dominant system are always larger
than the queues lengths in our system if both, the dominant system and our
system, start from the same initial state and have the same arrivals and
encounter the same packet collisions.

Figure 2: Queue 1, $Q_{1}$, Markov chain model for Dominant System 1.

For the first Dominant System, we assume that queue 2 will always have packets
to transmit; even if the queue was empty dummy packets will be transmitted
from queue 2. Clearly this will set a dominant system to our system since the
transmission of dummy packets can only result in more collisions and packet
losses. If for a given arrival rate pair ($\lambda_{1}$, $\lambda_{2}$) the
first dominant system is stable then clearly our system will be stable.
Therefore, the stability region of the first dominant system will provide an
inner bound for our system stability region.

For queue 1, the Markov chain describing the evolution of the queue is shown
in Fig. 2. Note that the Markov chain has two classes of states, namely,
$k_{F}$ and $k_{R}$ and $k=0,1,2,\cdots$. The subscript $F$ denotes first
transmission states and the subscript $R$ denotes retransmission states. Note
that in the retransmission states, queue 1 packet will always be delivered
since there is no collisions in these states (queue 2 is backing off); in
these states, either queue 1 length decreases by 1 if no arrival occurs or the
queue length will remain the same if an arrival occurs while being in these
retransmission states since the packet on the head of the queue is
successfully transmitted with probability 1.

The stability condition for queue 1 in Dominant System 1 is given in the
following lemma, which is proved in Appendix A.

###### Lemma III.1

The arrival rates for queue 1 and queue 2 in Dominant System 1 must satisfy
the following two conditions, respectively,

$\begin{split}\lambda_{1}&<\frac{p_{1}}{1+p_{1}p_{2}}\\\
\lambda_{2}&<p_{2}(1-\lambda_{1}-\lambda_{1}p_{2})\end{split}$ (1)

for the system to be stable.

###  III-B Dominant System 2

In the second Dominant System, we assume that queue 1 always has packets to
send (dummy packets are sent if the queue decides to transmit while being
empty). Again, this will decouple the interaction of the two queues since the
service rate of queue 2 will be independent of the state of queue 1.

The Markov chain for the evolution of queue 2 is shown in Fig. 3. Two classes
of states are defined in Fig. 3 and denoted by the subscripts ON and OFF. The
ON states denote the states where queue 2 can access the channel. The OFF
states denote the back off states where queue 1 is retransmitting its collided
packets. Note that the transitions from the $k_{\text{OFF}}$ state can be
either to the $k_{\text{ON}}$ state, if no arrival occurs in the slot, or to
the $(k+1)_{\text{ON}}$ state, if one arrival occurs in the slot. The OFF
states can never make a transition to a state with a lower number of packets
since in the OFF states queue 2 is in the back off mode and no access is
attempted.

The stability condition for queue 2 in Dominant System 2 is given in the
following lemma, which is proved in Appendix B (the analysis in Appendix B
will be based on the theory of homogeneous quasi birth-and-death (QBD) Markov
chains [9]).

###### Lemma III.2

The arrival rates for queue 1 and queue 2 in Dominant System 2 must satisfy
the following two conditions, respectively,

$\begin{split}\lambda_{1}&<\frac{p_{1}(1-p_{1}-\lambda_{2}p_{1})}{(1-p_{1})}\\\
\lambda_{2}&<\frac{p_{2}(1-p_{1})}{1+p_{1}p_{2}}\end{split}$ (2)

for the system to be stable.

Figure 3: Queue 2, $Q_{2}$, Markov chain model for Dominant System 2.

Note that the intersection of the two stability regions described in Lemma
III.1 and Lemma III.2 for a given access vector
$\mathbf{p}=[p_{1}\;p_{2}]^{T}$ (grey area in Fig. 4) can be interpreted as
follows. Define a new Dominant System (Dominant System 3) in which every queue
has always a packet to transmit. In this case, the transmission state of queue
1 can be represented by the two-state Markov chain model shown in Fig 5(a);
note that in this case queue 1 will be either in the “Transmission” state
denoted by $F$ or in the “Retransmission” state denoted by $R$ in Fig. 5(a).
Fig. 5(b) shows the Markov chain model for queue 2. Queue 2 will have two
states denoted by ON when queue 1 is in the F state and OFF when queue 1 is
the $R$ state (when queue 1 is in the $R$ state queue 2 will be in the back
off, OFF state). It is straightforward to show that the steady state
distributions for the two Markov chains shown in Fig. 5 are given by

$\begin{split}\pi_{F}&=\pi_{\text{ON}}=\frac{1}{1+p_{1}p_{2}}\\\
\pi_{R}&=\pi_{\text{OFF}}=\frac{p_{1}p_{2}}{1+p_{1}p_{2}}.\end{split}$ (3)

Figure 4: The union of the stability regions for the two dominant systems for
fixed access probabilities $p_{1}$ and $p_{2}$.

The service rate for queue 1 in Dominant System 3, $\mu_{1}^{\prime\prime}$,
is given by

$\begin{split}\mu_{1}^{\prime\prime}=p_{1}(1-p_{2})\pi_{F}+\pi_{R}=\frac{p_{1}}{1+p_{1}p_{2}},\end{split}$
(4)

where queue 1 is served with probability $p_{1}(1-p_{2})$ in the $F$ state and
with probability 1 in the $R$ state.

The service rate for queue 2 in Dominant System 3, $\mu_{2}^{\prime\prime}$,
is given by

$\begin{split}\mu_{2}^{\prime\prime}=p_{2}(1-p_{1})\pi_{\text{ON}}+0\times\pi_{\text{OFF}}=\frac{p_{2}(1-p_{1})}{1+p_{1}p_{2}},\end{split}$
(5)

where queue 2 is served with probability $p_{2}(1-p_{1})$ in the ON state and
with probability 0 in the OFF state.

(a) The two-state Markov chain model for queue 1 transmission state in
Dominant System 3.

(b) The two-state Markov chain model for queue 2 transmission state in
Dominant System 3.

Figure 5: Dominant System 3 Markov chain model.

###  III-C The Stability Region of the Random Access Protocol with Priorities

In this section, we derive the expression for the stability region of the
random access scheme with feedback exploitation where a priority is set to one
of the two queues. The following Lemma characterizes the stability region for
fixed random access probabilities, $p_{1}$ and $p_{2}$, for queue 1 and queue
2, respectively.

###### Lemma III.3

For a fixed random access probability vector $\mathbf{p}=[p_{1}\;p_{2}]^{T}$,
the stability region $\mathcal{R}(\mathbf{p})$ of the random access with
priorities is the union of the two regions described by

$\begin{split}\lambda_{2}<p_{2}(1-\lambda_{1}-\lambda_{1}p_{2})\end{split}$
(6)

when

$\begin{split}\lambda_{1}<\frac{p_{1}}{1+p_{1}p_{2}}\end{split}$ (7)

and

$\begin{split}\lambda_{1}<\frac{p_{1}(1-p_{1}-\lambda_{2}p_{1})}{(1-p_{1})}\end{split}$
(8)

when

$\begin{split}\lambda_{2}<\frac{p_{2}(1-p_{1})}{1+p_{1}p_{2}}.\end{split}$ (9)

for the system to be stable.

###### Proof:

The result in Lemma III.3 can be proved using the tool of stochastic dominance
presented in [6]. The indistinguishability argument at the stability region
boundary states that if the original system is unstable then its queues will
saturate and they will always have packets to transmit; therefore at the
boundaries of the stability region of the original system, the original system
will be indistinguishable from the dominant system and thus has the same
stability region boundaries [6]. ∎

The next theorem characterizes the entire stability region for the random
access protocol with priorities.

###### Theorem III.4

The boundary of the stability region, $\mathcal{R}$, of the random access
protocol with priorities, which is defined as the union of the
$\mathcal{R}(\mathbf{p})$ regions for the different
$\mathbf{p}=[p_{1}\;p_{2}]^{T}$ as

$\begin{split}\mathcal{R}=\bigcup_{\mathbf{p}\in[0,1]^{2}}\mathcal{R}(\mathbf{p})\end{split}$
(10)

can be characterized as

$\begin{split}\lambda_{2}=\left\\{\begin{array}[]{ll}1-2\lambda_{1}&\lambda_{1}\leq\frac{1}{3}\\\
\frac{(1-\lambda_{1})^{2}}{4\lambda_{1}}&\lambda_{1}>\frac{1}{3}.\end{array}\right.\end{split}$
(11)

###### Proof:

First, we will derive the boundary of the stability region defined in lemma
III.1, which can be found as

$\begin{split}\lambda_{2}^{*}(\lambda_{1})=&\mathbf{max}_{p_{1},p_{2}}\;p_{2}(1-\lambda_{1}-\lambda_{1}p_{2})\\\
&\text{{subject to }}0\leq p_{1}\leq 1,\;0\leq p_{2}\leq
1,\;\lambda_{1}<\frac{p_{1}}{1+p_{1}p_{2}}.\end{split}$ (12)

Ignoring the constraints in the last optimization problem and differentiating
the cost function in the last expression with respect to $p_{2}$ and equating
the derivative to 0 we can get the optimal value for $p_{2}$, denoted by
$p_{2}^{*}$, as222it is straightforward to prove that the cost function is
concave in $p_{2}$.

$p_{2}^{*}=\frac{1-\lambda_{1}}{2\lambda_{1}}.$ (13)

Note that for $\lambda_{1}\geq\frac{1}{3}$, we have $p_{2}^{*}\leq 1$. Also,
for $p_{1}=1$ and $p_{2}^{*}=\frac{1-\lambda_{1}}{2\lambda_{1}}$, the maximum
value for the first queue arrival rate is
$\frac{p_{1}}{1+p_{1}p_{2}}=\frac{2\lambda_{1}}{1+\lambda_{1}}>\lambda_{1}$
(i.e., the last constraint, $\lambda_{1}<\frac{p_{1}}{1+p_{1}p_{2}}$ is
satisfied with $p_{1}=1$), which means that for $\lambda_{1}\geq\frac{1}{3}$,
the value for $p_{2}$ that maximizes $\lambda_{2}$ for a given $\lambda_{1}$
is given by $p_{2}^{*}=\frac{1-\lambda_{1}}{2\lambda_{1}}$, with all the
constraints in (12) not being violated.

For $\lambda_{1}<\frac{1}{3}$, following similar steps to the
$\lambda_{1}\geq\frac{1}{3}$ case, we can easily prove that the value for
$p_{2}$ that maximizes $\lambda_{2}$ is giving by $p_{2}^{*}=1$; clearly the
values of $p_{1}=1$ and $p_{2}^{*}=1$ can be easily checked to satisfy the
constraints in (12) for $\lambda_{1}<\frac{1}{3}$.

Substituting the optimal values for $p_{2}$ for the different ranges of
$\lambda_{1}$ we can easily get the boundary of the stability region spanned
by the expression in lemma III.1 to be given by

$\begin{split}\lambda_{2}=\left\\{\begin{array}[]{ll}1-2\lambda_{1}&\lambda_{1}\leq\frac{1}{3}\\\
\frac{(1-\lambda_{1})^{2}}{4\lambda_{1}}&\lambda_{1}>\frac{1}{3}.\end{array}\right.\end{split}$
(14)

Finally, following a similar approach to that considered here it is
straightforward to show that the boundary derived in (14) is the boundary of
the stability regions defined in lemma III.2, which completes the proof. ∎

In Fig. 6, we have plotted the regions $\mathcal{R}(\mathbf{p})$, for $p_{1}$
and $p_{2}$ ranging from 0 to 1 with a step of 0.01, along with the derived
stability region boundary given in the previous theorem. Fig. 6 also shows the
stability region of the random access scheme, whose boundary is given by the
following relation [6]

$\sqrt{\lambda_{1}}+\sqrt{\lambda_{2}}=1.$ (15)

In Fig. 6, we also show the boundary of the stability region for the time
division (TD) based scheme (genie-aided), which serves as the stability region
upper bound, given by333Time Division (TD) corresponds to full coordination
between the two queues and requires knowledge of the queues arrival rates a
priori before dividing the resources (time slots in this case).

$\lambda_{1}+\lambda_{2}=1.$ (16)

Figure 6: The stability regions for the Random Access, Random Access with
Priorities, and Time Division schemes.

It is clear, and straightforward to analytically prove from the closed-form
stability region boundary expressions, that the stability region for the RA
scheme with priorities encloses the stability region of the RA scheme. This
can explained as follows. For a given arrival rate at the first queue,
$\lambda_{1}$, the RA with priority scheme will provide a better service rate
to that queue if compared to the RA scheme and this means that queue 1 will be
empty with a higher probability and this means that queue 2 will have a higher
service rate as well under the RA with priority scheme as compared to the RA
scheme. So setting a priority to the first queue in the retransmission will
also result in a service rate improvement for the second queue; this is
because the RA with priority scheme has some form of coordination between the
two queues in the retransmission stage. Allowing for collision free
retransmission from the first queue will decrease the amount of expected
collisions between the transmissions of the two queues and this will result in
better service rates for the two queues.

##  IV Conclusions

In this paper, we consider the problem of deriving the stability region for
random access protocol with feedback exploitation. We consider the case of two
interacting queue with priority set to one of the two queues. The two queues
will access the channel through a conventional random access protocol and in
the case of collision the higher priority queue will access the channel in the
next slot with probability 1 while the other queue will back off. We derive
the stability region for the random access with priorities protocol and prove
that it contains the stability region for the conventional random access
protocol. We show that not only the service rate for the higher priority queue
is enhanced but also the service rate for the other queue is improved if
compared to the conventional random access protocol.

##  Appendix A Proof of Lemma III.1

In this Appendix, we provide a proof for Lemma III.1. We start by calculating
the steady state distribution for the Markov chain shown in Fig. 2.

First, it is clear that $\epsilon_{0}=0$ since the queue can never be in a
retransmission state while being empty. Writing the balance equation around
$1_{R}$, we have

$\epsilon_{1}=\lambda_{1}p_{1}p_{2}\pi_{0}+\left(1-\lambda_{1}\right)p_{1}p_{2}\pi_{1}.$
(17)

Then around $0_{F}$, we have

$(\lambda_{1}p_{1}p_{2}+\lambda_{1}(1-p_{1}))\pi_{0}=\left(1-\lambda_{1}\right)\epsilon_{1}+\left(1-\lambda_{1}\right)p_{1}(1-p_{2})\pi_{1}.$
(18)

Substituting for $\epsilon_{1}$ from (17) into (18), and after some
manipulations, we can get

$\pi_{1}=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}\pi_{0}.$
(19)

Substituting from (19) into (17), we get

$\epsilon_{1}=\frac{\lambda_{1}p_{2}}{1-\lambda_{1}p_{2}}\pi_{0}.$ (20)

Writing the balance equation around $1_{F}$, we have

$\begin{split}&\left(1-\lambda_{1}p_{1}\left(1-p_{2}\right)-\left(1-\lambda_{1}\right)(1-p_{1})\right)\pi_{1}=\\\
&\quad\lambda_{1}\pi_{0}+\lambda_{1}\epsilon_{1}+\left(1-\lambda_{1}\right)\epsilon_{2}+\left(1-\lambda_{1}\right)p_{1}\left(1-p_{2}\right)\pi_{2}.\end{split}$
(21)

Around $2_{R}$, we have

$\epsilon_{2}=\lambda_{1}p_{1}p_{2}\pi_{1}+\left(1-\lambda_{1}\right)p_{1}p_{2}\pi_{2}.$
(22)

To get the relation between $\pi_{1}$ and $\pi_{2}$, we can substitute for the
values of $\epsilon_{1}$, $\pi_{0}$ and $\epsilon_{2}$ from equations (17),
(18) and (22), respectively in equation (21); after some tedious manipulation,
we get

$\pi_{2}=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}\pi_{1}.$
(23)

Substituting from (23) into (22), we get

$\epsilon_{2}=\frac{\lambda_{1}p_{2}}{1-\lambda_{1}p_{2}}\pi_{1}.$ (24)

Note that the Markov chain is repeating from stage 2 till the end. For $k\geq
2$, we have the following relations.

$\pi_{k}=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}\pi_{k-1}.$
(25) $\epsilon_{k}=\frac{\lambda_{1}p_{2}}{1-\lambda_{1}p_{2}}\pi_{k-1}.$ (26)

The last relation can be used to prove the following relation between
$\epsilon_{k}$ and $\epsilon_{k-1}$.

$\epsilon_{k}=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}\epsilon_{k-1}.$
(27)

The steady state distribution can now be written as follows.

  * •

$\epsilon_{0}=0$.

  * •

$\pi_{k}=\rho^{k}\pi_{0}$, $k\geq 1$ and
$\rho=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}$.

  * •

$\epsilon_{1}=\frac{\lambda_{1}p_{2}}{1-\lambda_{1}p_{2}}\pi_{0}$.

  * •

$\epsilon_{k}=\rho^{k-1}\epsilon_{1}$, $k\geq 2$.

This steady state distribution can be easily checked to satisfy the balance
equation at any general state (details are omitted since it is a rather
straightforward, yet very tedious, procedure).

To get the value of the steady state probabilities, we apply the following
normalization requirement.

$\begin{split}&\sum_{k=0}^{\infty}(\pi_{k}+\epsilon_{k})=1\\\
&\quad\quad\rightarrow\pi_{0}+\sum_{k=1}^{\infty}(\pi_{k}+\epsilon_{k})=\pi_{0}\left(1+\frac{\lambda_{1}p_{2}}{1-\lambda_{1}p_{2}}\right)\sum_{k=0}^{\infty}\rho^{k}=1,\end{split}$
(28)

where
$\rho=\frac{\lambda_{1}\left(1-p_{1}+\lambda_{1}p_{1}p_{2}\right)}{p_{1}\left(1-\lambda_{1}\right)\left(1-\lambda_{1}p_{2}\right)}$
as defined above.

Note that for the steady state distribution to exist, i.e. to have $\pi_{0}$
to be non zero, then we must have $\rho<1$, which is the stability condition
for queue 1 in this dominant system. Therefore, the stability condition can be
stated as

$\rho<1\;\rightarrow\lambda_{1}<\frac{p_{1}}{1+p_{1}p_{2}}.$ (29)

From the normalization condition in (28), we can get the value of $\pi_{0}$ as

$\pi_{0}=\frac{p_{1}-\lambda_{1}(1+p_{1}p_{2})}{p_{1}(1-\lambda_{1})}.$ (30)

In Dominant System 1, queue 2 will be served only in the states denoted by the
subscript $F$ in Fig. 2 since in the retransmission states, denoted by the
subscript $R$ in Fig. 2, queue 2 will be in the back off mode. Hence, the
service rate, $\mu_{2}$, for queue 2 in Dominant System 1 is given by

$\begin{split}\mu_{2}&=p_{2}(1-\lambda_{1})\pi_{0}+p_{2}(1-p_{1})\lambda_{1}\pi_{0}+\sum_{k=1}^{\infty}p_{2}(1-p_{1})\pi_{k}\\\
&=p_{2}(1-p_{1}\lambda_{1})\pi_{0}+\sum_{k=1}^{\infty}p_{2}(1-p_{1})\pi_{k},\end{split}$
(31)

where in the $0_{F}$ state, and with the arrival at the beginning of the slot
assumption, queue 2 is served with a rate of $p_{2}(1-\lambda_{1})\pi_{0}$
with no arrival at the beginning of the slot since queue 1 will not attempt
any random access since it is empty, and $p_{2}(1-p_{1})\lambda_{1}\pi_{0}$
with arrival at the slot beginning; for the other first transmission states,
queue 2 will be served if it decides to access the medium, which occurs with
probability $p_{2}$, and queue 1 decides not to access the medium, which
occurs with probability $(1-p_{1})$. After some manipulation, we can write the
expression for $\mu_{2}$ as

$\mu_{2}=p_{2}(1-\lambda_{1}-\lambda_{1}p_{2}).$ (32)

For the stability of queue 2, we must have

$\lambda_{2}<\mu_{2}=p_{2}(1-\lambda_{1}-\lambda_{1}p_{2}).$ (33)

##  Appendix B Proof of Lemma III.2

In this Appendix, we provide a proof for Lemma III.2. We start by calculating
the steady state distribution for the Markov chain shown in Fig. 3. The state
transition matrix, $\mathbf{\Phi}$, of the Markov chain shown in Fig. 3 can be
written as

$\mathbf{\Phi}=\left(\begin{array}[]{
ccccc}\mathbf{B}&\mathbf{A}_{0}&\mathbf{0}&\mathbf{0}&\cdots\\\
\mathbf{A}_{2}&\mathbf{A}_{1}&\mathbf{A}_{0}&\mathbf{0}&\cdots\\\
\mathbf{0}&\mathbf{A}_{2}&\mathbf{A}_{1}&\mathbf{A}_{0}&\cdots\\\
\mathbf{0}&\mathbf{0}&\mathbf{A}_{2}&\mathbf{A}_{1}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)$ (34)

where

$\begin{split}\mathbf{B}&=\left(\begin{array}[]{
cc}(1-\lambda_{2})+\lambda_{2}(1-p_{1})p_{2}&0\\\ 0&0\end{array}\right),\\\
\mathbf{A}_{0}&=\left(\begin{array}[]{ cc}(1-\lambda_{2})(1-p_{1})p_{2}&0\\\
0&0\end{array}\right),\\\ \mathbf{A}_{1}&=\left(\begin{array}[]{
cc}\lambda_{2}p_{2}(1-p_{1})+(1-\lambda_{2})(1-p_{2})&1-\lambda_{2}\\\
(1-\lambda_{2})p_{1}p_{2}&0\end{array}\right),\\\
\mathbf{A}_{2}&=\left(\begin{array}[]{ cc}\lambda_{2}&\lambda_{2}\\\
\lambda_{2}p_{1}p_{2}&\lambda_{2}\end{array}\right).\end{split}$

The steady state distribution vector is given by
$\mathbf{v}=[\pi_{0}^{\prime}\;\epsilon_{0}^{\prime}\;\pi_{1}^{\prime}\;\epsilon_{1}^{\prime}\;\pi_{2}^{\prime}\;\epsilon_{2}^{\prime}\;\cdots]^{T}$
and $\mathbf{v}=\mathbf{\Phi}\mathbf{v}$.

The state transition matrix $\mathbf{\Phi}$ is a block-tridiagonal matrix;
therefore the Markov chain shown in Fig. 3 is a homogeneous quasi birth-and-
death (QBD) Markov chain [9]. Note that to make the state transition matrix a
block-tridiagonal matrix we have added a transition from the $0_{\text{OFF}}$
state to the $1_{\text{ON}}$ state as shown in Fig. 7 and this will preserve
the structure of the state transitions between the different stages in the
Markov chain. Note that adding this transition will not affect the stationary
state distribution of the Markov chain as well as the balance equations since
$\epsilon_{0}^{\prime}=0$ even with the added transition since the Markov
chain will never enter the $0_{\text{OFF}}$ state444The analysis presented
here could have been used for analyzing the Markov chain shown in Fig. 2;
however, the structure of this Markov chain allowed for the use of a simpler
approach that was adopted in Appendix A.

Figure 7: The queue 2 Markov chain with added transition between
$0_{\text{OFF}}$ and $1_{\text{ON}}$ to make the state transition matrix a
block-tridiagonal matrix.

Define the vector
$\mathbf{v}_{k}^{\prime}=[\pi_{k}^{\prime}\;\epsilon_{k}^{\prime}]^{T}$. Note
that $\mathbf{v}_{0}^{\prime}=[\pi_{0}^{\prime}\;0]^{T}$. The steady state
distribution of the Markov chain shown in Fig. 3 satisfies the following
equation [9]

$\mathbf{v}^{\prime}_{k}=\mathbf{R}^{k}\mathbf{v}_{0}^{\prime},\quad k\geq 1,$
(35)

where the $2\times 2$ matrix $\mathbf{R}$ is given by the solution to the
following equation.

$\mathbf{A}_{0}+\mathbf{R}(\mathbf{A_{1}}-\mathbf{I}_{2})+\mathbf{R}^{2}\mathbf{A_{2}}=\mathbf{0}_{2\times
2},$ (36)

where $\mathbf{I}_{2}$ is the $2\times 2$ identity matrix and
$\mathbf{0}_{2\times 2}$ is the all zeros $2\times 2$ matrix.

To get the stationary distribution, we have to find the matrix

$\mathbf{R}=\left(\begin{array}[]{ cc}r_{11}&r_{12}\\\
r_{21}&r_{22}\end{array}\right).$

Note that for $\mathbf{v}_{1}^{\prime}=\mathbf{R}\mathbf{v}_{0}^{\prime}$,
where $\mathbf{v}_{0}^{\prime}=[\pi_{0}^{\prime}\;0]^{T}$ and
$\mathbf{v}_{1}^{\prime}=[\pi_{1}^{\prime}\;\epsilon_{1}^{\prime}]^{T}$.
Therefore, we have

$\begin{split}r_{11}=\frac{\pi_{1}^{\prime}}{\pi_{0}^{\prime}}\text{ and
}r_{21}=\frac{\epsilon_{1}^{\prime}}{\pi_{0}^{\prime}}.\end{split}$ (37)

Writing the balance equation around the $0_{\text{ON}}$ in Fig. 3, we have

$\begin{split}&(\lambda_{2}p_{1}p_{2}+\lambda_{2}(1-p_{2}))\pi_{0}^{\prime}=(1-\lambda_{2})(1-p_{1})p_{2}\pi_{1}^{\prime}\\\
&\quad\rightarrow\pi_{1}^{\prime}=\frac{\lambda_{2}(1-p_{2}+p_{1}p_{2})}{(1-\lambda_{2})(1-p_{1})p_{2}}\pi_{0}^{\prime}.\end{split}$
(38)

Therefore, we have

$\begin{split}r_{11}=\frac{\lambda_{2}(1-p_{2}+p_{1}p_{2})}{(1-\lambda_{2})(1-p_{1})p_{2}}.\end{split}$
(39)

Writing the balance equation around $1_{\text{OFF}}$, we have

$\epsilon_{1}^{\prime}=\lambda_{2}p_{1}p_{2}\pi_{0}^{\prime}+(1-\lambda_{2})p_{1}p_{2}\pi_{1}^{\prime}\rightarrow\epsilon_{1}^{\prime}=\frac{\lambda_{2}p_{1}}{1-p_{1}}\pi_{0}^{\prime}.$
(40)

Therefore, we have

$r_{21}=\frac{\lambda_{2}p_{1}}{1-p_{1}}.$ (41)

To get the values of $r_{12}$ and $r_{22}$, we consider the transition across
the border shown in Fig. 8. For the Markov chain to be positive recurrent then
the probability of going across the border in both directions must be the same
[10]; hence, we have

$\begin{split}(1-\lambda_{2})(1-p_{1})p_{2}\pi_{2}^{\prime}=(\lambda_{2}p_{1}p_{2}+\lambda_{2}(1-p_{2}))\pi_{1}^{\prime}+\lambda_{2}\epsilon_{1}^{\prime}.\end{split}$
(42)

But we have $\mathbf{v}_{2}^{\prime}=\mathbf{R}\mathbf{v}_{1}^{\prime}$, from
which we have
$\pi_{2}^{\prime}=r_{11}\pi_{1}^{\prime}+r_{12}\epsilon_{1}^{\prime}$; using
(39) and (42), we can easily show that

$\begin{split}r_{12}=\frac{\lambda_{2}}{(1-\lambda_{2})(1-p_{1})p_{2}}.\end{split}$
(43)

Figure 8: The segment of the Markov chain used to calculate the values of
$r_{12}$ and $r_{22}$.

Finally, the balance equation around $2_{\text{OFF}}$ can be written as

$\epsilon_{2}^{\prime}=\lambda_{2}p_{1}p_{2}\pi_{1}^{\prime}+(1-\lambda_{2})p_{1}p_{2}\pi_{2}^{\prime}=r_{21}\pi_{1}^{\prime}+r_{22}\epsilon_{2}^{\prime}.$
(41)

Substituting for $\pi_{2}^{\prime}$ from (42), we can easily show that

$r_{22}=\frac{\lambda_{2}p_{1}}{1-p_{1}}.$ (42)

Now the matrix $\mathbf{R}$ is given by

$\mathbf{R}=\left(\begin{array}[]{cc}\frac{\lambda_{2}(1-p_{2}+p_{1}p_{2})}{(1-\lambda_{2})(1-p_{1})p_{2}}&\frac{\lambda_{2}}{(1-\lambda_{2})(1-p_{1})p_{2}}\\\
\frac{\lambda_{2}p_{1}}{1-p_{1}}&\frac{\lambda_{2}p_{1}}{1-p_{1}}\end{array}\right),$
(43)

which can be easily checked to satisfy the balance equation given by (36).

To get the stationary distribution of the Markov chain shown in Fig. 3, we
apply the following normalization requirement.

$\sum_{k=0}^{\infty}(\pi_{k}^{\prime}+\epsilon_{k}^{\prime})=1\rightarrow[1\;1]\left(\sum_{k=0}^{\infty}\mathbf{R}^{k}\right)\mathbf{v}_{0}^{\prime}=1.$
(44)

For the summation $\left(\sum_{k=0}^{\infty}\mathbf{R}^{k}\right)$ to converge
we must have the spectral radius of the matrix $\mathbf{R}$,
$\textbf{sp}(\mathbf{R})$, to be less than one [9]555The spectral radius of a
matrix is the maximum over the magnitudes of its eigenvalues.. From (43), we
can easily get $\textbf{sp}(\mathbf{R})$ to be given by

$\begin{split}\textbf{sp}(\mathbf{R})=\frac{\lambda_{2}\left(1-p_{2}-\lambda_{2}p_{1}p_{2}+2p_{1}p_{2}+\sqrt{1-2p_{2}+p_{2}^{2}+4p_{1}p_{2}-2\lambda_{2}p_{1}p_{2}-2\lambda_{2}p_{1}p_{2}^{2}+\lambda_{2}^{2}p_{1}^{2}p_{2}^{2}}\right)}{2p_{2}\left(1-\lambda_{2}-p_{1}+\lambda_{2}p_{1}\right)}.\end{split}$
(45)

The requirement that $\textbf{sp}(\mathbf{R})<1$ can be used in the last
expression to get the stability condition of the second queue arrival rate
$\lambda_{2}$ as

$\begin{split}\lambda_{2}<\frac{p_{2}(1-p_{1})}{1+p_{1}p_{2}}.\end{split}$
(46)

Going back to the normalization requirement in (44), we have

$[1\;1]\left(\sum_{k=0}^{\infty}\mathbf{R}^{k}\right)\mathbf{v}_{0}^{\prime}=[1\;1]\left(\mathbf{I}_{2}-\mathbf{R}\right)^{-1}\left[\begin{array}[]{c}\pi_{0}^{\prime}\\\
0\end{array}\right]=1.$ (47)

From the last expression, we can easily prove that $\pi_{0}^{\prime}$ is given
by

$\pi_{0}^{\prime}=\frac{p_{2}-\lambda_{2}-p_{1}p_{2}-\lambda_{2}p_{1}p_{2}}{(1-\lambda_{2})(1-p_{1})p_{2}}.$
(48)

Note that the requirement that $\pi_{0}^{\prime}>0$, i.e. a non-zero empty
queue probability, is satisfied if
$\lambda_{2}<\frac{p_{2}(1-p_{1})}{1+p_{1}p_{2}}$, which is the queue
stability condition.

The service rate, $\mu_{1}^{\prime}$, for queue 1 in Dominant System 2 can now
be expressed as

$\begin{split}\mu_{1}^{\prime}&=p_{1}(1-\lambda_{2})\pi_{0}^{\prime}+p_{1}(1-p_{2})\lambda_{2}\pi_{0}^{\prime}+p_{1}(1-p_{2})\sum_{k=2}^{\infty}\pi_{k}^{\prime}+\sum_{k=2}^{\infty}\epsilon_{k}^{\prime}\\\
&=p_{1}(1-\lambda_{2}p_{2})\pi_{0}^{\prime}+[p_{1}(1-p_{2})\;1]\left(\sum_{k=1}^{\infty}\mathbf{R}^{k}\right)\left[\begin{array}[]{c}\pi_{0}^{\prime}\\\
0\end{array}\right]\\\
&=p_{1}(1-\lambda_{2}p_{2})\pi_{0}^{\prime}+[p_{1}(1-p_{2})\;1]\;\mathbf{R}\left(\mathbf{I}_{2}-\mathbf{R}\right)^{-1}\left[\begin{array}[]{c}\pi_{0}^{\prime}\\\
0\end{array}\right]\\\
&=\frac{p_{1}(1-p_{1}-\lambda_{2}p_{1})}{(1-p_{1})},\end{split}$ (49)

where in the OFF states, queue 1 is served with probability 1 since queue 2
will be in the back off mode. For the stability of queue 1 in Dominant System
2 we must have

$\lambda_{1}<\mu_{1}^{\prime}=\frac{p_{1}(1-p_{1}-\lambda_{2}p_{1})}{(1-p_{1})},$
(50)

which completes the proof.

## References

  * [1] A.K. Sadek, K.J.R. Liu, and Anthony Ephremides,  “Cognitive multiple access via cooperation: Protocol design and performance analysis,”  Information Theory, IEEE Transactions on, vol. 53, no. 10, pp. 3677–3696, 2007.
  * [2] S. Kompella, Gam D. Nguyen, J.E. Wieselthier, and Anthony Ephremides,  “Stable throughput tradeoffs in cognitive shared channels with cooperative relaying,”  in INFOCOM, 2011 Proceedings IEEE, 2011, pp. 1961–1969.
  * [3] Jeongho Jeon and Anthony Ephremides,  “The stability region of random multiple access under stochastic energy harvesting,”  in Information Theory Proceedings (ISIT), 2011 IEEE International Symposium on, 2011, pp. 1796–1800.
  * [4] A. Fanous and Anthony Ephremides,  “Effect of secondary nodes on the primary’s stable throughput in a cognitive wireless network,”  in Information Theory Proceedings (ISIT), 2012 IEEE International Symposium on, 2012, pp. 1807–1811.
  * [5] A. Fanous and A. Ephremides,  “Stable throughput in a cognitive wireless network,”  Selected Areas in Communications, IEEE Journal on, vol. 31, no. 3, pp. 523–533, 2013.
  * [6] R. Rao and A. Ephremides,  “On the Stability of Interacting Queues in a Multi-Access System,”  IEEE Trans. Info. Theory, vol. 34, pp. 918–930, Sept. 1988.
  * [7] W. Szpankowski,  “Stability Conditions for Some Multiqueue Distributed System: Buffered Random Access Systems,”  Adv. Appl. Probab., vol. 26, pp. 498–515, 1994.
  * [8] R. M. Loynes,  “The Stability of a Queue with Non-Independent Interarrival and Service Times,”  Proc. Cambridge Philos. Soc., pp. 497–520, 1962.
  * [9] G. Latouche and V. Ramaswami,  Introduction to Matrix Analytic Methods in Stochastic Modeling,  ASA-SIAM Series on Statistics and Applied Probability. Society for Industrial and Applied Mathematics, 1999.
  * [10] R.G. Gallager,  Discrete stochastic processes, vol. 101,  Kluwer Academic Publishers, 1996.


# Closest Wannier functions to a given set of localized orbitals

Taisuke Ozaki  Institute for Solid State Physics, The University of Tokyo,
Kashiwa 277-8581, Japan

###### Abstract

A non-iterative method is presented to calculate the closest Wannier functions
(CWFs) to a given set of localized guiding functions, such as atomic orbitals,
hybrid atomic orbitals, and molecular orbitals, based on minimization of a
distance measure function. It is shown that the minimization is directly
achieved by a polar decomposition of a projection matrix via singular value
decomposition, making iterative calculations and complications arising from
the choice of the gauge irrelevant. The disentanglement of bands is inherently
addressed by introducing a smoothly varying window function and a greater
number of Bloch functions, even for isolated bands. In addition to atomic and
hybrid atomic orbitals, we introduce embedded molecular orbitals in molecules
and bulks as the guiding functions, and demonstrate that the Wannier
interpolated bands accurately reproduce the targeted conventional bands of a
wide variety of systems including Si, Cu, the TTF-TCNQ molecular crystal, and
a topological insulator of Bi2Se3. We further show the usefulness of the
proposed method in calculating effective atomic charges. These numerical
results not only establish our proposed method as an efficient alternative for
calculating WFs, but also suggest that the concept of CWF can serve as a
foundation for developing novel methods to analyze electronic structures and
calculate physical properties.

##  I INTRODUCTION

The Wannier functions (WFs) Wannier37 ; Kohn59 ; Marzari12  play a central
role in analyzing electronic structures of real materials and advancing
electronic structure methods alongside density functional theory (DFT)
Hohenberg1964  and the other electronic structure theories Hamann09 ;
Lechermann06 . A widely adopted method for calculating WFs involves maximizing
their localization, which can be reformulated as minimizing the spread
function characterizing the variance of WFs in real space Marzari97 . The
concept of maximal localization leads to an elegant formulation, resulting in
maximally localized Wannier functions (MLWFs), and the compact representation
of the Hamiltonian that follows allows us to analyze the electronic structures
of real materials using a localized orbital picture. However, it is also
recognized that the minimization of the spread function can often encounter
local minima, particularly in large-scale systems having complicated
electronic structures Marzari12 . To overcome this challenge, methods aimed at
automated high-throughput Wannierisation have been developed Vitale20 ;
Damle18 ; Qiao23 . In addition, the other methods for generating
(nonorthogonal) WFs have been proposed based on projection methods Ku02 ;
Lu2004 ; Qian08 , which produce compact atomic like orbitals. The
orthogonalization of such nonorthogonal orbitals obtained by the projection
can be achieved by the Löwdin orthogonalization procedure Lowdin1950 ;
Cloizeaux1964 . These orthogonalized orbitals are utilized as an initial guess
of WFs in the minimization of the spread function to obtain MLWFs Marzari97 ;
Mostofi2008 . Apart from the role of initial guess, it is worth noting that
the Löwdin orthogonalized orbitals possess an important variational property,
which in particular minimizes the sum of squared distances between the
orthogonal and original nonorthogonal orbitals in the Hilbert space Carlson57
. As long as the original nonorthogonal orbitals are localized, the
localization of the Löwdin orthogonalized orbitals in real space is guaranteed
in a sense of minimization of the distance. In quantum chemistry, this
property is leveraged to develop methods for generating localized orthogonal
orbitals. These methods, which include the natural atomic and bond orbital
method Reed88  and the intrinsic atomic orbital method Knizia13 , involve
obtaining orthogonalized orbitals through Löwdin type orthogonalization,
giving heavy weight to occupied states via the density matrix. However, in
solid state physics the variational property inherent in the Löwdin
orthogonalization has not been fully explored yet. A reason for this is that
the Löwdin orthogonalization requires calculation of $S^{-1/2}$ for the
overlap matrix $S$, which causes numerical instability for nearly ill-
conditioned matrices. Due to the difficulty, advancements in methodologies
along these lines appear to be slow. Our study aims to develop an efficient
and robust method to generate the closest Wannier functions (CWFs) to a given
set of localized guiding functions. We will demonstrate that exploiting the
variational property leads to a versatile method that eliminates the need for
iterative optimization and avoids complications related to gauge choice.

The structure of this paper is as follows: In Section II, we present the
theory of calculating CWFs. Section III provides a detailed discussion on the
implementation of the method. A series of benchmark calculations are presented
in Section IV. Finally, in Section V, we summarize the theory of CWFs and
suggest a potential role for CWFs as a foundation for developing novel
methods.

##  II THEORY

Let us start by introducing Bloch functions $\\{\phi\\}$, which can be
obtained by solving the Kohn-Sham (KS) equation Kohn1965  within the density
functional theory (DFT) Hohenberg1964 , normalized as

$\displaystyle\langle\phi_{{\bf k}_{1}\mu}|\phi_{{\bf k}_{2}\nu}\rangle=N_{\rm
BC}\delta_{{\bf k}_{1}{\bf k}_{2}}\delta_{\mu\nu},$ (1)

where ${\bf k}$ and $\mu$ are the k-vector and the band index, respectively,
and $N_{\rm BC}$ is the number of primitive cells in the Born-von Karman (BvK)
boundary condition. We consider a projection of a localized guiding function
$Q$ onto the Bloch functions $\\{\phi\\}$ weighted with a window function
$w(\varepsilon)$ as

$\displaystyle|L_{{\bf R}p}\rangle=\frac{1}{N_{\rm BC}}\sum_{{{\bf
k},\mu}}{\rm e}^{-{\rm i}{\bf k}\cdot{\bf R}}|\phi_{{\bf k}\mu}\rangle a_{{\bf
k}\mu,p}$ (2)

with

$\displaystyle a_{{\bf k}\mu,p}=w(\varepsilon_{{\bf k}\mu})\langle\phi_{{\bf
k}\mu}|Q_{{\bf 0}p}\rangle,$ (3)

where ${\bf R}$ and $p$ are a translational lattice vector and the index of
localized orbital, respectively, $\varepsilon_{{\bf k}\mu}$ is the eigenvalue
of the KS equation, and $Q_{{\bf 0}p}\equiv Q_{{\bf R}p}~{}({\bf R}={\bf 0})$.
In the summation of Eq. (2), the numbers of k-point and the Bloch functions to
be included are $N_{\rm BC}$ and $N_{\rm band}$ per k-point, respectively.
Throughout the paper we do not consider the spin dependency on the formulation
for sake of simplicity, but the generalization is straightforward. The
localized functions $\\{Q_{{\bf R}p}\\}$, such as atomic orbitals, hybrid
orbitals, or molecular orbitals (MOs), are assumed to localize in the unit
cell specified by ${\bf R}$. When the window function $w(\varepsilon)$ is
taken to be unity, Eq. (2) is nothing but an expansion of $Q_{{\bf R}p}$ with
Bloch functions from the identity relation $\frac{1}{N_{\rm BC}}\sum_{{{\bf
k}\mu}}|\phi_{{\bf k}\mu}\rangle\langle\phi_{{\bf k}\mu}|=\hat{I}$. To
introduce an expansion of $Q_{{\bf R}p}$ by the Bloch functions $\\{\phi_{{\bf
k}\mu}\\}$ in a subspace, we use the following window function
$w(\varepsilon)$:

$\displaystyle
w(\varepsilon)=\frac{1-\exp(x_{0}+x_{1})}{\left(1+\exp(x_{0}))(1+\exp(x_{1})\right)}+\delta$
(4)

with the definition of $x_{0}\equiv\beta(\varepsilon_{0}-\varepsilon)$ and
$x_{1}\equiv\beta(\varepsilon-\varepsilon_{1})$, where $\beta=\frac{1}{k_{\rm
B}T}$ with the Boltzmann constant of $k_{\rm B}$ and a temperature of $T$. The
window function of Eq. (4) is obtained by subtracting 1 from the sum of two
functions of $1/(1+\exp(x))$. The last term $\delta$ is a small constant,
e.g., $10^{-12}$, which is introduced to avoid the ill-conditioning of the
matrix consisting of $a_{{\bf k}\mu,p}$ by Eq. (3) as discussed later on. As
shown in Fig. 1, the parameters of $\varepsilon_{0}$ and $\varepsilon_{1}$
($\varepsilon_{0}<\varepsilon_{1}$) determine the range of energy where the
Bloch states are included in the expansion of Eq. (2) with a large weight, and
the temperature $T$ gives the degree of smearing around $\varepsilon_{0}$ and
$\varepsilon_{1}$ in $w(\varepsilon)$. To control the degree of smearing, we
introduce the temperature of $T$ so that one can intuitively understand the
degree of smearing. It is noted that the choice of the window function by Eq.
(4) is not unique, and the other choices should give almost equivalent results
as long as they are chosen properly.

Figure 1:  Window functions of Eq. (4) with $k_{\rm B}T$=0.1, 1.0, and 3.0 eV.
$\varepsilon_{0}$ and $\varepsilon_{1}$ were set to be -5.0 and 5.0 eV,
respectively, and $\delta=10^{-12}$ was used.

The function $L_{{\bf R}p}$ defined by Eq. (2) must be similar to the
localized function $Q_{{\bf R}p}$, but they are non-orthogonal to each other
in general. We now consider generating a set of closest Wannier functions
(CWFs) $\\{W_{{\bf R}p}\\}$ to a set of the localized functions $\\{L_{{\bf
R}p}\\}$ in a sense that the sum of the squared distance between $L$ and $W$
in the Hilbert space is minimized. As well as Eq. (2), noting that the CWFs
can be expressed Marzari12  as

$\displaystyle|W_{{\bf R}p}\rangle=\frac{1}{N_{\rm BC}}\sum_{{{\bf
k},\mu}}{\rm e}^{-{\rm i}{\bf k}\cdot{\bf R}}|\phi_{{\bf k}\mu}\rangle b_{{\bf
k}\mu,p},$ (5)

and defining the residual function $R$:

$\displaystyle|R_{{\bf R}p}\rangle=|L_{{\bf R}p}\rangle-|W_{{\bf R}p}\rangle,$
(6)

the distance measure (DM) function $F$ we minimize is given by

$\displaystyle F[B]$ $\displaystyle=$ $\displaystyle\sum_{p}\langle R_{{\bf
0}p}|R_{{\bf 0}p}\rangle,$ (7) $\displaystyle=$ $\displaystyle\frac{1}{N_{\rm
BC}}\sum_{{\bf k}}X[B,{\bf k}]$

with

$\displaystyle X[B,{\bf k}]={\rm tr}\left[\left(A^{{\dagger}}({\bf
k})-B^{{\dagger}}({\bf k})\right)\left(A({\bf k})-B({\bf k})\right)\right],$
(8)

where the elements of the matrices $A({\bf k})$ and $B({\bf k})$ are given by
$a_{{\bf k}\mu,p}$ and $b_{{\bf k}\mu,p}$, respectively. The second line of
Eq. (7) is obtained by considering the orthonormality of Eq. (1), and
$X[B,{\bf k}]$ of Eq. (8) is regarded as the squared Frobenius norm of
$\left(A({\bf k})-B({\bf k})\right)$ Golub1996 . The same Bloch functions as
in Eq. (2) are included in the summation of Eq. (5). The number of CWFs to be
generated is $N_{\rm CWF}$ per unit cell, being equivalent to the number of
the guiding functions per cell, and $N_{\rm CWF}$ should be smaller than or
equal to $N_{\rm band}$ to guarantee the linear independency of the subspace
spanned by the CWFs, resulting in size of $N_{\rm band}\times N_{\rm CWF}$ for
the matrices $A({\bf k})$ and $B({\bf k})$. Assuming $B^{{\dagger}}({\bf
k})B({\bf k})=I$, where $I$ is of size $N_{\rm CWF}\times N_{\rm CWF}$, the
CWFs are ensured to form a set of orthonormal functions. Under this
constraint, we consider the optimization of the DM function $F$. Furthermore,
the matrix $B({\bf k})$ is regarded as a part of a unitary matrix of size
$N_{\rm band}\times N_{\rm band}$, and will be referred to as a partial
unitary matrix in subsequent discussions.

The minimum of the DM function of Eq. (7) is obtained by choosing $B({\bf k})$
as $U({\bf k})$ which is calculated by the polar decomposition of $A({\bf
k})=U({\bf k})P({\bf k})$, where $U({\bf k})$ and $P({\bf k})$ are unitary and
hermitian, respectively Fan1955 . Let us prove the statement below. The polar
decomposition of $A({\bf k})$ is obtained via the singular value decomposition
(SVD) of $A({\bf k})$ as

$\displaystyle A=W\Sigma V^{{\dagger}}=WV^{{\dagger}}V\Sigma
V^{{\dagger}}=UP,$ (9)

where we dropped the dependency on ${\bf k}$ for simplicity of the notation,
and hereafter we will denote the dependency if necessary. $W$ and $V$ are the
left and right singular matrices in size of $N_{\rm band}\times N_{\rm CWF}$
and $N_{\rm CWF}\times N_{\rm CWF}$, respectively, and $\Sigma$ is the
singular value diagonal matrix in size of $N_{\rm CWF}\times N_{\rm CWF}$.
Note that $U\equiv WV^{{\dagger}}$ and $P\equiv V\Sigma V^{{\dagger}}$. It is
worth mentioning that $W$ and $U$ are partial unitary matrices, and hold
$W^{{\dagger}}W=I$ ($WW^{{\dagger}}\neq I$) and $U^{{\dagger}}U=I$
($UU^{{\dagger}}\neq I$), and that $V$ is a full unitary matrix holding
$V^{{\dagger}}V=I$ ($VV^{{\dagger}}=I$). We evaluate $X$ of Eq. (8) for each
${\bf k}$ with both the matrix $U$ obtained by the polar decomposition and an
arbitrary partial unitary matrix $B$, and calculate the difference as

$\displaystyle X[U]-X[B]$ $\displaystyle=$ $\displaystyle 2{\rm
tr}\left[\frac{1}{2}\left(A^{{\dagger}}B+B^{{\dagger}}A\right)-P\right],$ (10)
$\displaystyle=$ $\displaystyle 2{\rm tr}\left[\Sigma D-\Sigma\right],$
$\displaystyle=$ $\displaystyle 2\sum_{n}\sigma_{n}\left(d_{nn}-1\right)$

with

$\displaystyle
D=\frac{1}{2}\left(V^{{\dagger}}U^{{\dagger}}BV+V^{{\dagger}}B^{{\dagger}}UV\right).$
(11)

The second line of Eq. (10) is derived by using the polar decomposition of $A$
and $P=V\Sigma V^{{\dagger}}$, and $\sigma_{n}$ and $d_{nn}$ in the third line
are diagonal elements of $\Sigma$ and $D$, respectively. The upper bound of
the diagonal elements of the hermitian matrix $D$ is found to be unity, since
the matrix $D$ consists of the sum of the product of partial unitary matrices
of $V^{{\dagger}}U^{{\dagger}}$ ($V^{{\dagger}}B^{{\dagger}}$) and $BV$
($UV$). Another proof for the upper bound is also given based on the Cauchy-
Schwarz inequality in the appendix. By considering the upper bound of the
diagonal elements of $D$ and $0\leq\sigma_{n}$, the third line of Eq. (10)
leads to the following inequality:

$\displaystyle X[U]\leq X[B].$ (12)

The equality of Eq. (12) holds if $B=U$. If some of $\sigma_{n}$ is zero, the
corresponding $d_{nn}$ in Eq. (10) can be arbitrarily chosen. So, we see from
Eq. (11) that $U$ is not uniquely determined. If all the singular values of
$\sigma_{n}$ are positive, all the $d_{nn}$s should be unity when $X[U]=X[B]$.
The uniqueness of the solution can be confirmed as follows: Since the matrix
$D$ consists of the products of the partial unitary matrices as discussed
above, the case that all the diagonal elements $d_{nn}$ are unity happens when

$\displaystyle U^{{\dagger}}B+B^{{\dagger}}U=2I,$ (13)

which is derived by equating $D=I$ and multiplying $V$ and $V^{{\dagger}}$
from the left and right of the equation, respectively. By noting again that
$B$ and $U$ in Eq. (13) are the partial unitary matrices, it is found that Eq.
(13) holds if and only if $B=U$. Thus, we have proven the statement that the
minimum of the DM function of Eq. (7) is uniquely determined when $B=U$ as
long as all the singular values of $\sigma_{n}$ are positive, since $F[B]$ of
Eq. (7) consists of the sum of $X[B,{\bf k}]$ over k. The uniqueness of $U$
itself is related to that of the SVD for the matrix $A$ of Eq. (3). If $A$ has
$N_{\rm CWF}$ positive singular values which are non-degenerate, the SVD is
uniquely determined except for the non-unique phases in $W$ and $V$. Even if
there are degenerate singular values, the matrix $U$ is invariant, since the
freedom of the unitary transformation $K$ for the degenerate left and right
singular vectors is canceled out as
$U=WKK^{{\dagger}}V^{{\dagger}}=WV^{{\dagger}}$ with $KK^{{\dagger}}=I$.
Further, noting that the non-unique phases in $W$ and $V$ are canceled out
when the matrix $U$ is computed as $WV^{{\dagger}}$, we conclude that $U$ is
uniquely determined if $A$ has the $N_{\rm CWF}$ positive singular values. The
violation from the positive definiteness of the singular values can be avoided
by the small constant of $\delta$ in Eq. (4). On the other hand, if $\delta$
in Eq. (4) is set to be 0 and a narrow window, including fewer than $N_{\rm
CWF}$ eigenstates for a given ${\bf k}$, is used with a small $k_{\rm B}T$,
the matrix $A$ has $N_{\rm PSV}$ positive singular values, where $N_{\rm
PSV}<N_{\rm CWF}$. In this case, a subspace of the $\left(N_{\rm CWF}-N_{\rm
PSV}\right)$ dimension is arbitrary chosen to form $W$ and $V$ of the
dimension of $N_{\rm CWF}$ in addition to the subspace of the $N_{\rm PSV}$
dimension, which breaks the symmetry of CWFs. The situation can be avoided by
introducing a small constant $\delta$ in the window function of Eq. (4). The
simple treatment restores the symmetry of CWFs, allowing us to calculate
symmetry preserving CWFs for a wide variety of choices in the energy range of
the window function. It is also noted that the non-unique phase of the Bloch
functions is canceled out via the polar decomposition of $A$ and Eq. (5).
Thus, we see that the proposed method is free from complications arising from
the choice of gauge.

Using Eqs. (8) and (9), the minimum of the DM function $F$ at $B=U$ is given
by

$\displaystyle F[U]=\frac{1}{N_{\rm BC}}\sum_{{\bf k},p}\left(\sigma_{{\bf
k}p}-1\right)^{2},$ (14)

where $\sigma_{{\bf k}p}$ is a singular value of the matrix $A({\bf k})$. From
Eq. (14), we find that the mean squared distance between $L$ and $W$ is
related to the deviation of $\sigma$ from unity.

The proposed method based on the polar decomposition of $A$ is closely related
to the Löwdin orthogonalization Lowdin1950  as shown below. The Fourier
transform of the overlap integrals for $\\{L\\}$ is given by

$\displaystyle\sum_{{\bf R}}{\rm e}^{{\rm i}{\bf k}\cdot{\bf R}}\langle
L_{{\bf 0}p}|L_{{\bf R}q}\rangle=\sum_{\mu}a_{{\bf k}\mu,p}^{*}a_{{\bf
k}\mu,q},$ (15)

which is obtained by noting that $\frac{1}{N_{\rm BC}}\sum_{{\bf R}}{\rm
e}^{{\rm i}({\bf k}-{\bf k}^{\prime})\cdot{\bf R}}=\delta_{{\bf k}{\bf
k}^{\prime}}$. By writing Eq. (15) as $S({\bf k})=A^{{\dagger}}({\bf k})A({\bf
k})$ in a matrix form, and using Eq. (9), we have the following relation:

$\displaystyle S({\bf k})=A^{{\dagger}}({\bf k})A({\bf k})=V^{{\dagger}}({\bf
k})\Sigma^{2}({\bf k})V({\bf k}).$ (16)

Comparing Eq. (16) with $P({\bf k})=V({\bf k})\Sigma({\bf
k})V^{{\dagger}}({\bf k})$, one obtains a relation $P({\bf k})=S^{1/2}({\bf
k})$. If $P({\bf k})$ is invertible, the matrix $U({\bf k})$ is given with Eq.
(9) by

$\displaystyle U({\bf k})=A({\bf k})S^{-1/2}({\bf k}).$ (17)

The matrix $U$ obtained by Eq. (17) is exactly equivalent to that by the
Löwdin orthogonalization Lowdin1950 , and the closest property of the Löwdin
orthogonalized orbitals to a given set of orbitals was proven in Ref.
Carlson57 . On the other hand, the proposed method does not require the
calculation of $S^{-1/2}$, and can be applied in a numerical stable manner
even to the case that the matrix $A$ is nearly ill-conditioned. The definition
of $A$ by Eq. (3) with the window function allows us to calculate CWFs in
various respects, e.g., heavily weighting to occupied states and disentangling
of localized states embedded in the wider bands as demonstrated in Section IV.
In this sense, the method we propose can be regarded as a generalization of
the Löwdin orthogonalization.

The disentanglement of bands can be easily performed by properly selecting the
window function of Eq. (3). Since the $N_{\rm band}(\geq N_{\rm CWF})$ Bloch
functions per k-point are included in Eqs. (2) and (5), it should be
emphasized that the proposed method always disentangles $N_{\rm band}$ states
to generate $N_{\rm CWF}$ Wannier functions. A couple of examples will be
shown in Section IV, including valence and low-lying conduction bands for the
diamond Si and narrow $3d$-bands embedded in the wider $4s$-band for the face
centered cubic (FCC) Cu.

Once the $k$-dependent $U({\bf k})$ are obtained by the polar decomposition,
tight-binding (TB) parameters are calculated as expectation values of the KS
Hamiltonian $\hat{H}_{\rm KS}$ by summing contributions over ${\bf k}$ and
$\mu$ as

$\displaystyle t_{{\bf 0}p,{\bf R}q}$ $\displaystyle=$ $\displaystyle\langle
W_{{\bf 0}p}|\hat{H}_{\rm KS}|W_{{\bf R}q}\rangle,$ (18) $\displaystyle=$
$\displaystyle\sum_{{\bf k},\mu}\varepsilon_{{\bf k}\mu}u_{{\bf
k}\mu,p}^{*}u_{{\bf k}\mu,q}{\rm e}^{-{\rm i}{\bf k}\cdot{\bf R}},$

where $u_{{\bf k}\mu,q}$ is the matrix element of $U({\bf k})$. Like with
MLWFs, the TB parameters can be used for the Wannier interpolation in
calculations of band structures and physical quantities.

The computational procedure to generate the CWFs is summarized as follows:

  1. 1.

Determining $\varepsilon_{0}$, $\varepsilon_{1}$, and $k_{\rm B}T$ in Eq. (4)
by checking the band dispersion of a system of interest. The parameters should
be chosen properly so that the targeted eigenstates can be included.

  2. 2.

Choosing a set of localized orbitals {$Q$}. Atomic orbitals, hybrid orbitals,
and MOs might be possible choices. Depending on the symmetry of the targeted
eigenstates, a proper set of localized orbitals needs to be chosen, e.g.,
$d$-orbitals need to be employed in the generation of CWFs for the $d$-bands
in FCC Cu as shown later on.

  3. 3.

Calculation of $A({\bf k})$ by Eq. (3). The calculation of the overlap
$\langle\phi_{{\bf k}\mu}|Q_{{\bf 0}p}\rangle$ depends on the implementation
of the KS method.

  4. 4.

Performing the SVD of $A({\bf k})$. A proper numerical library might be used.

  5. 5.

Calculation of $U({\bf k})$. Using the left and right singular matrices,
$W({\bf k})$ and $V({\bf k})$, of $A({\bf k})$, $U({\bf k})$ is calculated as
$W({\bf k})V^{{\dagger}}({\bf k})$.

  6. 6.

Calculations of the CWFs and TB parameters. The CWFs can be calculated using
Eq. (5) with $B=U$ obtained by the step 5. The calculation is expressed by the
sum of the matrix-matrix product over ${\bf k}$. Also, the TB parameters are
obtained by Eq. (18).

It should be stressed that the minimization of $F$ can be efficiently
performed in a non-iterative manner by the steps 1 to 6 using the polar
decomposition via the SVD. The computational cost of each step is estimated to
be $O(N_{\rm BC}N_{\rm band}N_{\rm CWF}N_{\rm basis})$, $O(N_{\rm BC}N_{\rm
band}^{3})$, $O(N_{\rm BC}N_{\rm band}N_{\rm CWF}^{2})$, or $O(N_{\rm
BC}N_{\rm band}N_{\rm CWF}N_{\rm basis})$ for the step 3, 4, 5, or 6,
respectively, where $N_{\rm basis}$ is the number of basis functions to expand
the KS orbitals. If the localized basis functions are used, the computational
cost of the step 3 becomes $O(N_{\rm BC}N_{\rm band}N_{\rm CWF})$.

For isolated systems, the Brillouin zone sampling is limited to only the
$\Gamma$ point. No modification of the proposed method is needed to calculate
the CWFs.

The theoretical framework we have discussed so far is general for any
implementation of the DFT-KS method. However, the choice of the the localized
functions $\\{Q\\}$ in Eq. (3) may depend on the implementation. So, we will
discuss how the localized functions $\\{Q\\}$ can be properly generated in
Section III.

##  III IMPLEMENTATION

###  III.1 General

We have implemented the proposed method into the OpenMX DFT software package
OpenMX ; Ozaki2005 ; Duy2014  which is based on norm-conserving
pseudopotentials (PPs) MBK1993 ; Theurich2001  and optimized pseudo-atomic
orbitals (PAOs) Ozaki2003 ; Ozaki2004  as basis set. All the benchmark
calculations were performed with a computational condition of a production
level. The basis functions used are listed in Table 1. In the abbreviation of
basis functions such as H7.0-s2p2d1, H stands for the atomic symbol, 7.0 the
cutoff radius (Bohr) in the generation by the confinement scheme Ozaki2003 ;
Ozaki2004 , and s2p2d1 means the employment of two, two, and one optimized
radial functions for the $s$-, $p$-, and $d$-orbitals, respectively. The
radial functions were optimized by a variational optimization method Ozaki2003
. These basis functions we used can be regarded as at least double zeta plus
double polarization basis sets if we follow the terminology of Gaussian or
Slater-type basis functions. Valence states included in the PPs are listed in
Table 1. All the PPs and PAOs we used in the study were taken from the
database (2019) in the OpenMX website OpenMX , which were benchmarked by the
delta gauge method Lejaeghere2016 . Real space grid techniques were used for
the numerical integrations and the solution of the Poisson equation using fast
Fourier transform (FFT) with the energy cutoff of 250 to 500 Ryd Ozaki2005 .
We used a generalized gradient approximation (GGA) proposed by Perdew, Burke,
and Ernzerhof to the exchange-correlation functional Perdew1996 . An
electronic temperature of 300 K was used to count the number of electrons by
the Fermi-Dirac function for all the systems we considered. For all the
calculations, $\delta=10^{-12}$ was used in the window function of Eq. (4),
and $N_{\rm band}$ was set to be equilavent to the number of basis functions
in both the summation of Eqs. (2) and (5).

Table 1:  Basis functions and valence states included in PPs. ∗For the calculations of atomic charges, the basis functions listed in Table 3 were used.  Element |  |  | Basis functions |  |  | Valence states in PP
---|---|---|---|---|---|---
H |  |  | H7.0-s2p2d1 |  |  | $1s$
C |  |  | C6.0-s3p2d2 |  |  |  $2s$, $2p$
N |  |  | N6.0-s3p2d2 |  |  |  $2s$, $2p$
Na |  |  |  ∗Na9.0-s3p3d2f2 |  |  |  $2s$, $2p$, $3s$
Si |  |  | Si7.0-s2p2d2f1 |  |  |  $3s$, $3p$
Cl |  |  |  ∗Cl9.0-s3p3d2f2 |  |  |  $3s$, $3p$
S |  |  | S7.0-s3p2d2f1 |  |  |  $3s$, $3p$
Cu |  |  | Cu6.0-s3p3d3f1 |  |  |  $3s$, $3p$, $3d$, $4s$
Se |  |  | Se7.0-s3p2d2 |  |  |  $4s$, $4p$
Bi |  |  | Bi8.0-s3p2d2f1 |  |  |  $5d$, $6s$, $6p$

###  III.2 Choice of Guiding Functions $\\{Q\\}$

Depending on the targeted eigenstates in the window function of Eq. (4), one
can choose either atomic orbitals, hybrid atomic orbitals, or MOs as the
guiding functions $\\{Q\\}$ in Eq. (3). In this subsection we discuss the
three kinds of the localized guiding functions $\\{Q\\}$, i.e., atomic
orbitals, hybrid atomic orbitals, and embedded MOs in molecules and bulks, and
how they can be generated in our implementation.

In our implementation, the Bloch function $\phi_{{\bf k}\mu}$ is expanded by
PAOs $\chi$ Ozaki2003 ; Ozaki2004  as

$\displaystyle|\phi_{{\bf k}\mu}\rangle=\sum_{{\bf R}}{\rm e}^{{\rm i}{\bf
k}\cdot{\bf R}}\sum_{i\alpha}c_{{\bf k}\mu,i\alpha}|\chi_{{\bf
R}i\alpha}\rangle,$ (19)

where $i$ and $\alpha$ are atomic and orbital indices, and $c$ is a linear
combination of PAO (LCPAO) coefficient. Also note that $\langle{\bf
r}|\chi_{{\bf R}i\alpha}\rangle\equiv\chi_{i\alpha}({\bf r}-{\bf\tau}_{i}-{\bf
R})$, where ${\bf\tau}_{i}$ is the position of atom $i$. We use PAOs as the
guiding functions $\\{Q\\}$ of atomic orbitals, corresponding to valence
orbitals or specific orbitals such as localized $d$-orbitals in a narrow
energy window.

The choice of the atomic orbitals gives a good guiding function, however, the
resultant CWFs may not recover the symmetry of CWFs respecting the bond
direction to the neighboring atoms unlike hybrid atomic orbitals due to the
closest property of CWFs. In this case, it would be better to use the hybrid
atomic orbitals rather than the atomic orbitals as explained below. Let us
introduce a projection operator for the occupied space by

$\displaystyle\hat{P}=\frac{1}{N_{\rm BC}}\sum_{{\bf k}\mu}|\phi_{\bf
k\mu}\rangle f(\varepsilon_{\bf k\mu})\langle\phi_{\bf k\mu}|,$ (20)

where $f$ is the Fermi-Dirac function. With the projection operator, the
density matrix is calculated by

$\displaystyle\rho_{{\bf R}i\alpha,{\bf R}^{\prime}j\beta}$ $\displaystyle=$
$\displaystyle\sum_{{\bf k}\mu}\langle\widetilde{\chi}_{{\bf
R}i\alpha}|\hat{P}|\widetilde{\chi}_{{\bf R}^{\prime}j\beta}\rangle,$
$\displaystyle=$ $\displaystyle\frac{1}{N_{\rm BC}}\sum_{{\bf k}\mu}{\rm
e}^{{\rm i}{\bf k}\cdot\left({\bf R}-{\bf
R}^{\prime}\right)}f(\varepsilon_{\bf k\mu})c_{{\bf k}\mu,i\alpha}c^{*}_{{\bf
k}\mu,j\beta},$

where $\widetilde{\chi}$ is the dual orbital defined by

$\displaystyle|\widetilde{\chi}_{{\bf R}i\alpha}\rangle$ $\displaystyle=$
$\displaystyle\sum_{{\bf R}^{\prime}j\beta}|\chi_{{\bf
R}^{\prime}j\beta}\rangle S_{{\bf R}^{\prime}j\beta,{\bf R}i\alpha}^{-1}$ (22)

holding the following orthonormality:

$\displaystyle\langle\chi_{{\bf R}i\alpha}|\widetilde{\chi}_{{\bf
R}^{\prime}j\beta}\rangle$ $\displaystyle=$
$\displaystyle\langle\widetilde{\chi}_{{\bf R}i\alpha}|\chi_{{\bf
R}^{\prime}j\beta}\rangle=\delta_{{\bf
RR}^{\prime}}\delta_{ij}\delta_{\alpha\beta}.$ (23)

By setting ${\bf R}={\bf R}^{\prime}={\bf 0}$ and $i=j$ in Eq. (LABEL:eq:DM),
and diagonalizing the diagonal block element consisting of $\rho_{{\bf
0}i\alpha,{\bf 0}i\beta}$, which is associated with selected atomic orbitals
on the atomic site $i$, we obtain the hybrid atomic orbitals respecting the
symmetry of the bond direction to the neighboring atoms. We employ the hybrid
atomic orbitals as the localized functions $\\{Q\\}$ in our implementation.
When the same atomic orbitals are chosen to form the diagonal block element as
for the case of atomic orbitals $\\{Q\\}$, the resultant CWFs span the same
subspace in the Hilbert space.

For some systems the electronic structures can be rather easily understood by
employing MOs as building blocks. A molecular crystal is such a case, where
the eigenstates near the chemical potential can be approximately constructed
by a linear combination of the MOs of constituting molecules. Another example
is the bond in molecules and bulks. The bond between two atoms embedded in the
system might be analyzed by MOs associated with the two atoms. Here we show
how embedded MOs in molecules and bulks can be calculated in a simple
procedure. Let us start by noting that using Eq. (20) and assuming that the
Bloch function is expressed by Eq. (19) the total number of electrons for a
spin degenerate case is given by Ozaki2018

$\displaystyle N_{\rm ele}$ $\displaystyle=$ $\displaystyle 2{\rm
tr}\left[\hat{P}\right],$ (24) $\displaystyle=$ $\displaystyle 2\sum_{{\bf
R}i\alpha}\langle\widetilde{\chi}_{{\bf R}i\alpha}|\hat{P}|\chi_{{\bf
R}i\alpha}\rangle,$

In the second line of Eq. (24) we used the following identity operator:

$\displaystyle\hat{I}=\sum_{{\bf R}i\alpha}|\chi_{{\bf
R}i\alpha}\rangle\langle\widetilde{\chi}_{{\bf R}i\alpha}|=\sum_{{\bf
R}i\alpha}|\widetilde{\chi}_{{\bf R}i\alpha}\rangle\langle\chi_{{\bf
R}i\alpha}|.$ (25)

Since $N_{\rm ele}=N_{\rm BC}N_{\rm ele}^{(0)}$ with $N_{\rm
ele}^{(0)}=2\sum_{i\alpha}\langle\widetilde{\chi}_{{\bf
0}i\alpha}|\hat{P}|\chi_{{\bf 0}i\alpha}\rangle$, it is enough to consider
$N_{\rm ele}^{(0)}$. We utilize the formula for $N_{\rm ele}^{(0)}$ to
calculate embedded MOs in molecules and bulks, and rewrite it the sum of
partial traces as

$\displaystyle N_{\rm ele}^{(0)}=2\sum_{g}{\rm
tr}_{g}\left[(\widetilde{\chi}_{{\bf 0}g}|\hat{P}|\chi_{{\bf
0}g})\right]=2\sum_{g}{\rm tr}_{g}\left[\Lambda_{g}\right],$ (26)

where $g$ is the index of group including PAOs on grouped atoms, e.g., a set
of PAOs on a molecule. The notation of $|\chi_{{\bf 0}g})$ stands for a subset
of PAOs as

$\displaystyle|\chi_{{\bf 0}g})=\left(\cdots,|\chi_{{\bf
0}i1}\rangle,|\chi_{{\bf 0}i2}\rangle,\cdots,\right),$ (27)

where PAOs on all the atoms in the group $g$ are included in the subset. The
notation of $|\widetilde{\chi}_{{\bf 0}g})$ is the counterpart for the dual
orbitals. Using the identity operator of Eq. (25), $\Lambda_{g}$ in Eq. (26)
is given by

$\displaystyle\Lambda_{g}=\sum_{{\bf R}g^{\prime}}\rho_{{\bf 0}g,{\bf
R}g^{\prime}}S_{{\bf R}g^{\prime},{\bf 0}g}$ (28)

with definition of block elements:

$\displaystyle\rho_{{\bf R}g,{\bf R}^{\prime}g^{\prime}}$ $\displaystyle=$
$\displaystyle(\widetilde{\chi}_{{\bf R}g}|\hat{P}|\widetilde{\chi}_{{\bf
R}^{\prime}g^{\prime}}),$ (29) $\displaystyle S_{{\bf R}g,{\bf
R}^{\prime}g^{\prime}}$ $\displaystyle=$ $\displaystyle(\chi_{{\bf
R}g}|\chi_{{\bf R}^{\prime}g^{\prime}}).$ (30)

We now introduce the embedded MOs orbitals in molecules and bulks by a right-
singular vectors of $\Lambda_{g}$. Since the elements of $\Lambda_{g}$ are
real in case of the collinear DFT, the right-singular vectors can be obtained
with real components. However, the SVD of $\Lambda_{g}$ tends to produce the
right-singular vectors with complex components, which results in the complex
CWFs making analysis complicated. To obtain the real right-singular vectors,
we perform the eigendecomposition of $\Lambda_{g}^{{\dagger}}\Lambda_{g}$ as

$\displaystyle\Lambda_{g}^{{\dagger}}\Lambda_{g}=Y_{g}\Omega_{g}^{2}Y_{g}^{{\dagger}},$
(31)

where $\Omega_{g}^{2}$ is the diagonal matrix consisting of the eigenvalues of
$\Lambda_{g}^{{\dagger}}\Lambda_{g}$, and $Y_{g}$ is a unitary matrix
consisting of the corresponding eigenvectors $\\{y_{g,\nu}\\}$. By noting that
the SVD of $\Lambda_{g}$ is given by $Z_{g}\Omega_{g}Y_{g}^{{\dagger}}$, we
see that the right-singular vectors $Y_{g}$ can be obtained by diagonalizing
$\Lambda_{g}^{{\dagger}}\Lambda_{g}$ as $Y_{g}$. The square roots of the
eigenvalues of $\Lambda_{g}^{{\dagger}}\Lambda_{g}$ are singular values of
$\Lambda_{g}$, and can be approximately regarded as the occupation of the
corresponding eigenvector $y_{g,\nu}$. If $\Lambda_{g}$ is a real matrix,
$Y_{g}$ obtained by the diagonalization of
$\Lambda_{g}^{{\dagger}}\Lambda_{g}$ is guaranteed to be a real unitary
matrix. We further normalize the eigenvector $y_{g,\nu}$ by considering the
overlap integrals of Eq. (30) as

$\displaystyle|\bar{y}_{g,\nu}\rangle=\frac{|y_{g,\nu}\rangle}{\sqrt{\langle
y_{g,\nu}|S_{{\bf 0}g,{\bf 0}g}|y_{g,\nu}\rangle}},$ (32)

and calculate an expectation value of the KS Hamiltonian $\hat{H}_{\rm KS}$
with $\bar{y}_{g,\nu}$ as

$\displaystyle\epsilon_{g,\nu}=\langle\bar{y}_{g,\nu}|\hat{H}_{\rm
KS}|\bar{y}_{g,\nu}\rangle.$ (33)

After $\\{\bar{y}_{g,\nu}\\}$ are ordered based on the expectation values, we
employ selected ones among $\\{\bar{y}_{g,\nu}\\}$, e.g., ones near the
chemical potential, as the guiding functions $\\{Q\\}$ in Eq. (3) by
monitoring the expectation values and the corresponding singular values. They
are what we call embedded MOs in molecules and bulks. As demonstrated in
Section IV, $\\{\bar{y}\\}$ work as the good guiding functions to capture the
electronic structure of a molecular crystal.

##  IV NUMERICAL RESULTS

###  IV.1 Wannier Interpolated Bands

Figure 2:  Interpolated bands (red circles) of silicon in the diamond
structure, calculated by the TB Hamiltonian derived from CWFs with (a) $k_{\rm
B}T$=0.01 and (b) 3.0 eV in the window function of Eq. (4), respectively.
$\varepsilon_{0}$ and $\varepsilon_{1}$ relative to the chemical potential
were set to be -15.0 and 0.0 eV, respectively. The solid line is the original
one directly calculated by the conventional scheme. The number of k-points for
the Brillouin zone sampling was 13 $\times$ 13 $\times$ 13\. The experimental
lattice constant of 5.43 Å was used. The values of the DM function per CWF are
0.571 and 0.444 for (a) and (b), respectively.  Figure 3:  CWFs for (a) Si,
(b) Cu, (c) TTF in TTF-TCNQ, and (d) TCNQ in TTF-TCNQ. In all the cases,
isovalues of $\pm 0.04$ (orange:0.04, blue:-0.04) are used for drawing the
isosurfaces using OpenMX Viewer Lee2019 . The computational conditions for
(a), (b), (c), and (d) are the same as those in Fig. 2 (b), Fig. 4 (a), Fig. 5
(b), and Fig. 5 (b), respectively.

We present four numerical examples to demonstrate the broad applicability of
the proposed method across various systems.

In Figs. 2 (a) and (b), the interpolated bands of Si in the diamond structure,
calculated by the TB Hamiltonian derived from CWFs, are compared with those by
the conventional calculation. The hybrid atomic orbitals consisting of valence
minimal orbitals, which are one $s$\- and a set of $p_{x}$-, $p_{y}$-, and
$p_{z}$-orbitals per Si, were used as the guiding functions $\\{Q\\}$. For the
window function of Eq. (4) we use -15.0 and 0.0, relative to the chemical
potential, in eV for $\varepsilon_{0}$ and $\varepsilon_{1}$, respectively,
which covers the energy range of valence bands. On the other hand, two $k_{\rm
B}T$ of 0.01 and 3.0 eV were used to check how conduction bands are reproduced
depending on the effect of smearing. As clearly seen, the two cases well
reproduce the valence bands, and the conduction bands are also reproduced well
in the case of $k_{\rm B}T=3.0$ eV. Evaluating the DM function per CWF using
Eq. (14), we obtain values of 0.571 and 0.444 for $k_{\rm B}T$ of 0.01 and 3.0
eV, respectively. For $k_{\rm B}T=0.01$ eV, the first four singular values at
each $k$ point are close to 1, while the remaining four approach
$\delta(=10^{-12})$. Conversely, for $k_{\rm B}T=3.0$ eV, the eight singular
values at each $k$ point distribute ranging from 1.2 to 0.1. These differences
in the distribution of singular values explain why the value of the DM
function for $k_{\rm B}T=0.01$ eV is larger than that for $k_{\rm B}T=3.0$ eV.
One of the obtained CWFs is shown in Fig. 3 (a), which represents a $p$-like
CWF. In case of $k_{\rm B}T=0.01$ eV, the conduction bands are treated less
importantly due to heavily weighting to the valence bands, resulting in the
poor reproduction of the conduction bands. At first glance, one may consider
that this is an undesirable feature in such a treatment. However, the feature
enables us to generate the minimal atomic-like orthogonal orbitals, which well
span the subspace by the valence bands, and to calculate an effective charge
allocated to each atom using the minimal atomic-like orthogonal orbitals in an
unbiased way. We will demonstrate the calculation of the effective atomic
charges and stress the usefulness of the method later on.

Figure 4:  Interpolated bands (red circles) of copper in the FCC structure,
calculated by the TB Hamiltonian derived from CWFs with (a)
$\varepsilon_{0}=-5.5$, $\varepsilon_{1}=-1.0$ eV, $k_{\rm B}T=1.0$ eV, and
hybrid $3d$-orbitals as $\\{Q\\}$, and (b) $\varepsilon_{0}=-11.0$,
$\varepsilon_{1}=24.0$ eV, $k_{\rm B}T=3.0$ eV, and hybrid $3d$-, $4s$-,
$4p$-orbitals as $\\{Q\\}$. The solid line is the original one directly
calculated by the conventional scheme. The number of k-points for the
Brillouin zone sampling was 21 $\times$ 21 $\times$ 21\. The experimental
lattice constant of 3.61 Å was used. The values of the DM function per CWF are
0.081 and 0.211 for (a) and (b), respectively.

The disentanglement of bands in metals is demonstrated for Cu in the FCC
structure as shown in Figs. 4 (a) and (b). By selecting the parameters
$\varepsilon_{0}$ and $\varepsilon_{1}$, so that the $3d$ bands are included,
and only five $d$-orbitals as the guiding functions $\\{Q\\}$, the five
$d$-bands are reproduced as shown in Fig. 4 (a). We see that the $3d$-bands
are properly disentangled from the dispersive $4s$ band, and one of the
obtained CWFs preserves the shape of $d_{z^{2}}$-orbital as shown in Fig. 3
(b). When $3d$-, $4s$-, and $4p$-orbitals are used as the guiding functions
$\\{Q\\}$, a broad range of bands are reproduced as shown in Fig. 4 (b). Since
our PP of Cu includes $3s$, $3p$, $3d$, and $4s$ states, the original band
structure has the deep $3s$\- and $3p$-bands. We set the parameters in the
window function so as to discard the $3s$\- and $3p$-bands, and use the $3d$-,
$4s$-, and $4p$-orbitals, not $3s$\- and $3p$-orbitals, as the guiding
functions $\\{Q\\}$. So, the $3s$\- and $3p$-bands are not included in the
interpolated bands (not shown). The examples show that the CWFs can be
flexibly and easily constructed in accordance with the purpose of study
without numerical difficulties.

Table 2:  Expectation values $\epsilon$ (eV) of the KS Hamiltonian calculated by Eq. (33), which is relative to the chemical potential, and singular values $\omega$ of $\Lambda_{g}$ for embedded MOs in the TTF-TCNQ molecular crystal.  | TTF |  |  |  |  | TCNQ |
---|---|---|---|---|---|---|---
index | $\epsilon$ | $\omega$ |  |  | index | $\epsilon$ | $\omega$
23 | -0.069 | 1.359 |  |  | 34 | 0.135 | 1.082
24 | 0.719 | 1.346 |  |  | 35 | 0.375 | 1.329
25 | 0.761 | 1.449 |  |  | 36 | 0.477 | 1.336
26 | 3.408 | 0.864 |  |  | 37 | 4.928 | 0.535
27 | 5.947 | 0.047 |  |  | 38 | 9.389 | 0.036
28 | 6.197 | 0.047 |  |  | 39 | 9.795 | 0.019
Figure 5:  Interpolated bands, red lines in (a) and red circles in (b), of
TTF-TCNQ, calculated by the TB Hamiltonian derived from CWFs with (a)
$\varepsilon_{0}=-22.0$ and $\varepsilon_{1}=1.0$ eV relative to the chemical
potential, $k_{\rm B}T=2.0$ eV, and hybrid atomic valence orbitals as
$\\{Q\\}$, and (b) $\varepsilon_{0}=-1.0$ and $\varepsilon_{1}=1.0$ eV
relative to the chemical potential, $k_{\rm B}T=0.01$ eV, and the 26th MO and
the 37th MO for TTF and TCNQ as $\\{Q\\}$. The solid line is the original one
directly calculated by the conventional scheme. The number of k-points for the
Brillouin zone sampling was 5 $\times$ 19 $\times$ 3\. An experimental crystal
structure was used Kistenmacher1974 . The values of the DM function per CWF
are 0.441 and 0.022 for (a) and (b), respectively.

The interpolated bands for a molecular crystal, consisting of
tetrathiafulvalene (TTF) and tetracyanoquinodimethane (TCNQ), are shown in
Figs. 5 (a) and (b). By employing hybrid atomic orbitals as the guiding
functions $\\{Q\\}$, the wide range of bands are reproduced as shown in Fig. 5
(a). On the other hand, as shown in Fig. 5 (b), only four bands near the
chemical potential are reproduced using the embedded MOs in the bulk as the
guiding functions $\\{Q\\}$. The MOs were calculated by grouping the TTF and
TCNQ molecules separately as explained in Sec. III B. The expectation values
calculated by Eq. (33) and the singular values $\omega$ of $\Lambda_{g}$ are
listed in Table 2. It can be seen that the singular values largely change from
$\omega_{26}$ to $\omega_{27}$ and from $\omega_{37}$ to $\omega_{38}$ for TTF
and TCNQ, respectively. So, the 26th MO and the 37th MO for TTF and TCNQ,
respectively, were used as the guiding functions $\\{Q\\}$. Since there are
two TTF molecules and two TCNQ molecules in the unit cell, we have the four
embedded MOs as $\\{Q\\}$, resulting in the four bands. In Figs. 3 (c) and (d)
two of the resultant CWFs are shown, which localize in TTF and TCNQ molecules,
respectively, as expected. The value of the DM function per CWF is found to be
0.022, which implies that the guiding MOs are very close to the resultant
CWFs. The example demonstrates that the proposed method provides a direct way
to calculate CWFs for targeted bands in molecular crystals.

We extend the proposed method to the non-collinear DFT vonBarth1972 ;
Kubler1988  with spin-orbit interaction (SOI) Theurich2001 . The KS orbitals
are expressed by two-component spinor, and the SOI is introduced by fully
relativistic $j$-dependent PPs OpenMX . The theoretical framework is readily
extended to the non-collinear DFT with the SOI without any difficulty.
However, it should be noted that the hybrid atomic orbitals and embedded MOs
in molecules and bulks as the guiding functions $\\{Q\\}$ can be complex
functions in the case. In Figs. 6 (a) and (b), we show the Wannier
interpolated bands of Bi2Se3 which is known to be a topological insulator when
the SOI is included in the DFT calculation. For the cases without and with the
SOI, corresponding to Figs. 6 (a) and (b), respectively, it is confirmed that
the interpolated bands accurately reproduce the original ones regardless of
existence of the band inversion, demonstrating a wide variety of applicability
of the proposed method.

Figure 6:  Interpolated bands (red circles) of Bi2Se3 (a) without and (b) with
the spin-orbit interaction, calculated by the TB Hamiltonian derived from
CWFs. The solid line is the original ones directly calculated by the
conventional scheme. $\varepsilon_{0}=-6.0$ and $\varepsilon_{1}=3.0$ eV
relative to the chemical potential, $k_{\rm B}T=1.0$ eV, and hybrid Bi
$6p$-orbitals and Se $4p$-orbitals as $\\{Q\\}$, were used in both the
calculations. The number of k-points for the Brillouin zone sampling was 7
$\times$ 7 $\times$ 7\. An experimental crystal structure was used
Nakajima1963 . The values of the DM function per CWF are 0.163 and 0.184 for
(a) and (b), respectively.

###  IV.2 Effective Atomic Charges

From the construction, one can calculate CWFs closest to (hybrid) valence
atomic orbitals, while fully respecting the occupied states. The CWFs provide
us a unique way to evaluate effective atomic charges. Since the CWFs are a set
of orthonormal functions, the effective atomic charge $\kappa_{i}$ of atom $i$
is easily calculated without arbitrariness as

$\displaystyle\kappa_{i}$ $\displaystyle=$ $\displaystyle
N_{i}^{\rm(val)}-\sum_{p\in{\rm atom}~{}i}\langle W_{{\bf 0}p}|\hat{P}|W_{{\bf
0}p}\rangle,$ (34) $\displaystyle=$ $\displaystyle
N_{i}^{\rm(val)}-\frac{1}{N_{\rm BC}}\sum_{p\in{\rm atom}~{}i}\sum_{\bf
k,\mu}f(\varepsilon_{\bf k\mu})u_{{\bf k}\mu,p}^{*}u_{{\bf k}\mu,p},$

where $\hat{P}$ is the projection operator defined by Eq. (20), and
$N_{i}^{\rm(val)}$ is the number of valence electrons in the PP of atom $i$.
Though the summation over orbitals $p$ belonging to atom $i$ is taken into
account in Eq. (34), one can analyze the orbital resolved charges as well.
Table 3 shows effective atomic charges in a HCN molecule and the NaCl bulk,
calculated by Eq. (34) and the Mulliken population analysis Mulliken1955 .
Both the systems are known to be notorious cases because of the difficulty in
calculating the effective charges Knizia13 . We see that the effective charges
by the CWFs quickly converge as a function of basis functions, while those by
the Mulliken population are highly dependent on the choice of the basis
function. The effective atomic charges obtained from the CWFs could serve as a
valuable tool for analyzing electronic structures in a manner related to the
approach discussed in Ref. Knizia13 . Further investigation in this direction
will be conducted in future studies.

Table 3:  Effective atomic charges in a HCN molecule and the NaCl bulk, calculated by the proposed method (CWF) and Mulliken population analysis (MP) with a variety of basis functions, where $A$ represents the constituting atoms. Hybrid atomic valence orbitals are used as the guiding functions $\\{Q\\}$. For the window function, $\varepsilon_{0}=-55.0$ and $\varepsilon_{1}=0.0$ eV for HCN and $\varepsilon_{0}=-35.0$ and $\varepsilon_{1}=0.0$ eV for NaCl, relative to the chemical potential, and $k_{\rm B}T=0.1$ eV are used. The number of k-points for the Brillouin zone sampling is 9 $\times$ 9 $\times$ 9 for the NaCl bulk with the lattice constant of 5.63 Å.  HCN |  | CWF |  |  |  |  | MP |
---|---|---|---|---|---|---|---|---
Basis | H | C | N |  |  | H | C | N
$A$6.0-s1p1 | 0.077 | -0.052 | -0.025 |  |  | 0.384 | -0.164 | -0.221
$A$6.0-s2p2 | 0.069 | 0.003 | -0.073 |  |  | -0.070 | 0.321 | -0.252
$A$6.0-s2p2d1 | 0.066 | 0.009 | -0.076 |  |  | -0.008 | 0.393 | -0.385
$A$6.0-s3p3d2 | 0.066 | 0.008 | -0.074 |  |  | 0.110 | 0.298 | -0.408
$A$6.0-s3p3d2f2 | 0.066 | 0.008 | -0.074 |  |  | 0.167 | 0.297 | -0.464
NaCl | CWF |  |  |  |  | MP |  |
Basis | Na | Cl |  |  |  | Na | Cl |
$A$9.0-s2p1 | 0.391 | -0.391 |  |  |  | 0.716 | -0.716 |
$A$9.0-s3p2 | 0.422 | -0.422 |  |  |  | 0.595 | -0.595 |
$A$9.0-s3p3d2 | 0.421 | -0.421 |  |  |  | 0.158 | -0.158 |
$A$9.0-s3p3d2f2 | 0.421 | -0.421 |  |  |  | -0.062 | 0.062 |

##  V CONCLUSIONS

We presented a non-iterative method to calculate the closest Wannier functions
(CWFs) to a given set of localized guiding functions such as atomic orbitals,
hybrid atomic orbitals, and molecular orbitals (MOs). We defined the distance
measure (DM) function $F$ by the sum of the squared distance between the
projection functions $L$ and Wannier functions $W$ in the Hilbert space, and
considered minimizing the DM function. It was shown that the minimization of
the DM function is achieved by the polar decomposition of the projection
matrix $A$ with a window function via the singular value decomposition (SVD)
in a non-iterative manner. The CWFs can be uniquely constructed as long as the
projection matrix $A$ is positive definite. It was also shown that the method
is free from subtle choice of the gauge. The disentanglement of bands is
naturally taken into account by introducing a smoothly varying window
function, and including $N_{\rm band}$ Bloch functions to generate $N_{\rm
CWF}$ CWFs, where $N_{\rm CWF}\leq N_{\rm band}$. Even for isolated bands, the
method always performs the disentanglement of bands without any difficulty. We
have implemented the proposed method into the OpenMX code, which is based on
the density functional theory (DFT), numerical pseudo-atomic orbitals (PAOs),
and norm-conserving pseudopotentials (PPs), and introduced three types of
guiding functions, i.e., atomic orbitals, hybrid atomic orbitals, and embedded
MOs in molecules and bulks. The first two are easily employed from the PAOs.
For the last one, we developed a method to calculate embedded MOs in molecules
and bulks, which focuses on a partial trace formula of the projection operator
for the occupied space and applies SVD to the partial matrix for the
projection operator. The interpolated bands by tight-binding (TB) models
derived from the CWFs reproduce well the targeted conventional bands of a wide
variety of systems including Si, Cu, the TTF-TCNQ molecular crystal, and a
topological insulator of Bi2Se3. These successful reproduction of targeted
bands clearly demonstrates a wide variety of applicability of the proposed
method. We further show the usefulness of the proposed method in calculating
effective atomic charges, implying that the CWFs closest to atomic orbitals
can be used as a measure to analyze electronic structures from one system to
the others. Thus, we conclude that the proposed method is an alternative way
in efficiently calculating WFs, and the concept of CWFs will provide a basis
for development of novel methods of analyzing electronic structures and
calculating physical properties.

###### Acknowledgements.

The author wishes to express gratitude to Prof. Yoshiaki Sugimoto for
inspiring the development of CWFs through collaborative research. Thanks are
also extended to Mr. Ryotaro Koshoji for his insightful comments on the
theoretical aspect. Part of the computation in this study was carried out
using the computational facility of the Institute for Solid State Physics at
the University of Tokyo.

##  Appendix A The upper bound of the diagonal elements of the matrix $D$

A proof for the upper bound of the diagonal elements of the matrix $D$ is
given in the appendix. Note that $UV$ and $BV$ in Eq. (11) are partial unitary
matrices in size of $N_{\rm band}\times N_{\rm CWF}$. Since $D$ is hermitian,
its eigenvalues are real. Let $\lambda$ be an eigenvalue of $D$ and $x$ be the
corresponding eigenvector. Then, the following equation holds:

$\displaystyle\frac{1}{2}\left(V^{{\dagger}}U^{{\dagger}}BV+V^{{\dagger}}B^{{\dagger}}UV\right)|x\rangle=\lambda|x\rangle.$
(35)

Defining $|y\rangle=UV|x\rangle$ and $|z\rangle=BV|x\rangle$, and operating
$\langle x|$ from the left side of Eq. (35), we have

$\displaystyle\langle y|z\rangle+\langle z|y\rangle=2\lambda\langle
x|x\rangle.$ (36)

Noting $\langle y|z\rangle+\langle z|y\rangle\leq 2|\langle y|z\rangle|$, and
using the Cauchy-Schwarz inequality $|\langle y|z\rangle|\leq\sqrt{\langle
y|y\rangle\langle z|z\rangle}$, we obtain

$\displaystyle\langle y|z\rangle+\langle z|y\rangle\leq 2\sqrt{\langle
y|y\rangle\langle z|z\rangle}.$ (37)

Since $\langle y|y\rangle=\langle z|z\rangle=\langle x|x\rangle$, combining
Eq. (36) and Eq. (37) results in

$\displaystyle\lambda\langle x|x\rangle\leq\langle x|x\rangle.$ (38)

Considering $\langle x|x\rangle\neq 0$, the upper bound of the eigenvalues of
$D$ is found to be

$\displaystyle\lambda\leq 1.$ (39)

Noting that the matrix $D$ can be written by $x$ and $\lambda$ as

$\displaystyle D=\sum_{\nu}|x_{\nu}\rangle\lambda_{\nu}\langle x_{\nu}|,$ (40)

the diagonal elements $d_{nn}$ of the matrix $D$ is given by

$\displaystyle d_{nn}=\sum_{\nu}|\langle n|x_{\nu}\rangle|^{2}\lambda_{\nu},$
(41)

where $\langle n|x_{\nu}\rangle$ is the $n$-th element in the vector
$x_{\nu}$. Since the upper bound of $\lambda$ is unity, and $\\{x\\}$ forms a
unitary matrix, we obtain the upper bound of the diagonal elements of the
matrix $D$ as

$\displaystyle d_{nn}\leq 1.$ (42)

## References

  * (1) G.H. Wannier, Phys. Rev. 52, 191 (1937).
  * (2) W. Kohn, Phys. Rev. 115, 809 (1959).
  * (3) N. Marzari, A.A. Mostofi, J.R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * (4) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * (5) D.R. Hamann and D. Vanderbilt, Phys. Rev. B 79, 045109 (2009).
  * (6) F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. Posternak, A. Yamasaki, and O.K. Andersen, Phys. Rev. B 74, 125120 (2006).
  * (7) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * (8) V. Vitale, G. Pizzi, A. Marrazzo, J.R. Yates, N. Marzari, and A.A. Mostofi, npj Computational Materials 6, 66 (2020).
  * (9) A. Damle and L. Lin, Multiscale Model. Simul. 16, 1392 (2018).
  * (10) J. Qiao, G. Pizzi, and N. Marzari, arXiv:2306.00678.
  * (11) W. Ku, H. Rosner, W.E. Pickett, and R.T. Scalettar, Phys. Rev. Lett. 89, 167204 (2002).
  * (12) W.C. Lu, C.Z. Wang, T.L. Chan, K. Ruedenberg, and K.M. Ho, Phys. Rev. B 70, 041101(R) (2004).
  * (13) X. Qian, J. Li, L. Qi, C.-Z. Wang, T.-L. Chan, Y.-X. Yao, K.-M. Ho, and S. Yip, Phys. Rev. B 78, 245112 (2008).
  * (14) P.-O. Löwdin, J. Chem. Phys. 18, 365 (1950).
  * (15) J. Des Cloizeaux, Phys. Rev. 135, A685 (1964).
  * (16) A.A. Mostofi, J.R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
  * (17) B.C. Carlson and J.M. Keller, Phys. Rev. 105, 102 (1957).
  * (18) A.E. Reed, L.A. Curtiss, and F. Weinhold, Chem. Rev. 88, 899 (1988).
  * (19) G. Knizia, J. Chem. Theory Comput. 9, 4834 (2013).
  * (20) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * (21) Matrix Computations (3rd ed.). Baltimore: The Johns Hopkins University Press (1996), G. Golub and C.F Van Loan.
  * (22) K. Fan and A. J. Hoffman, Proc. Amer. Math. Soc. 6, 111 (1955).
  * (23) The code, OpenMX, pseudo-atomic basis functions, and pseudopotentials are available under terms of the GNU-GPL on a web site (http://www.openmx-square.org/).
  * (24) T. Ozaki and H. Kino, Phys. Rev. B 72, 045121 (2005).
  * (25) T.V.T. Duy and T. Ozaki, Comput. Phys. Commun. 185, 777 (2014).
  * (26) I. Morrison, D.M. Bylander, L. Kleinman, Phys. Rev. B 47, 6728 (1993).
  * (27) G. Theurich and N.A. Hill, Phys. Rev. B 64, 073106 (2001).
  * (28) T. Ozaki, Phys. Rev. B. 67, 155108, (2003).
  * (29) T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004).
  * (30) K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I.E. Castelli, S.J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J.K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J.A. Flores-Livas, K.F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E.K. Gross, A. Gulans, F. Gygi, D.R. Hamann, P.J. Hasnip, N.A. Holzwarth, D. Iu̧san, D.B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçü kbenli, Y.O. Kvashnin, I.L. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordströ m, T. Ozaki, L. Paulatto, C.J. Pickard, W. Poelmans, M.I. Probert, K. Refson, M. Richter, G.M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunströ m, A. Tkatchenko, M. Torrent, D. Vanderbilt, M.J. van Setten, V. Van Speybroeck, J.M. Wills, J.R. Yates, G.X. Zhang, and S. Cottenier, Science 351, aad3000 (2016).
  * (31) J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (32) T. Ozaki, M. Fukuda, and G. Jiang, Phys. Rev. B 98, 245137 (2018).
  * (33) Y.-T. Lee and T. Ozaki, J. Mol. Graph. Model. 89, 192 (2019); https://www.openmx-square.org/viewer/
  * (34) T.J. Kistenmacher, T.E. Phillips, and D.O. Cowan, Acta Cryst. B30, 763 (1974).
  * (35) S. Nakajima, J. Phys. Chem. Solids 24, 479 (1963).
  * (36) U. von Barth and L. Hedin, J. Phys. C: Solid State Phys. 5, 1629 (1972).
  * (37) J. Kubler, K.-H. Hock, J. Sticht, and A. R. Williams, J. Phys. F: Met. Phys. 18, 469 (1988).
  * (38) R.S. Mulliken, J. Chem. Phys. 23, 1833 (1955).